Generate Triangle Files for Medical Images Using VTK

The tcl script march.tcl uses the Visualization Toolkit to generate triangle files from cross-section data using the Marching Cubes algorithm. Make Your Own Visible Woman describes how to use this script. The script is broken into several sections:

Verify Parameters

Setup The Pipeline

Read The Volume Data

Generate Triangles

Decimate The Triangles

Smooth The Triangle Vertices

Generate Normals

Generate Triangles Strips

Write The Triangles

Execute The Pipeline

Set The Start Methods

Execute Each Pipeline Filter

This paper describes the script section by section.

Verify Parameters

Certain parameters are required and certain parameters have defaults. This section makes sure the required NAME, RESOLUTION, STUDY and END_SLICE parameters are defined. If a required parameter is missing, the script exits with an error message.

set OK 1
# Set defaults for any parameters not defined
if {[info exists NAME] == 0} {
    puts "march.tcl: NAME is not defined"
    set OK 0
}
if {[info exists RESOLUTION] == 0} {
    puts "march.tcl: RESOLUTION is not defined"
    set OK 0
}
if {[info exists STUDY] == 0} {
    puts "march.tcl: STUDY is not defined"
    set OK 0
}
if {[info exists END_SLICE] == 0} {
    puts "march.tcl: END_SLICE is not defined"
    set OK 0
}
if {[info exists VALUE] == 0} {
    puts "march.tcl: VALUE is not defined"
    set OK 0
}
if { $OK == 0 } {
    puts "march.tcl: One or more required parameters are missing!"
    exit
}

if {[info exists SLICE_ORDER] == 0} {set SLICE_ORDER si}
if {[info exists HEADER_SIZE] == 0} {set HEADER_SIZE 0}
if {[info exists PIXEL_SIZE] == 0} {set PIXEL_SIZE 1}
if {[info exists START_SLICE] == 0} {set START_SLICE 1}
if {[info exists BYTE_ORDER] == 0} {set BYTE_ORDER BigEndian}
if {[info exists REDUCTION] == 0} {set REDUCTION 1}
if {[info exists FEATURE_ANGLE] == 0} {set FEATURE_ANGLE 60}
if {[info exists DECIMATE_ANGLE] == 0} {set DECIMATE_ANGLE $FEATURE_ANGLE}
if {[info exists SMOOTH_ANGLE] == 0} {set SMOOTH_ANGLE $DECIMATE_ANGLE}
if {[info exists DECIMATE_ITERATIONS] == 0} {set DECIMATE_ITERATIONS 1}
if {[info exists SMOOTH_ITERATIONS] == 0} {set SMOOTH_ITERATIONS 10}
if {[info exists SMOOTH_FACTOR] == 0} {set SMOOTH_FACTOR .01}
if {[info exists DECIMATE_REDUCTION] == 0} {set DECIMATE_REDUCTION 1}
if {[info exists DECIMATE_ERROR] == 0} {set DECIMATE_ERROR .0002}
if {[info exists DECIMATE_ERROR_INCREMENT] == 0} {set DECIMATE_ERROR_INCREMENT .0002}
if {[info exists VOI] == 0} {
    set VOI "0 [expr $RESOLUTION - 1] 0 [expr $RESOLUTION - 1] $START_SLICE $END_SLICE]"
}
if {[info exists SAMPLE_RATE] == 0} {set SAMPLE_RATE "1 1 1"}

Setup The Pipeline

vtk uses a demand driven visualization pipeline to process data. Of course, you know that because you've read the book. The pipeline is discussed in Chapter 4. This section of the script defines the components of the pipeline to read the slices, create triangles, reduce their number, smooth the vertices and write the triangles to a file.

Read The Volume Data

We assume here that the all the data to be processed was acquired with a constant center landmark. In vtk, the origin of the data applies to the lower left of an image volume. We calculate the vtk origin such that the x,y center of the volume will be (0,0). We specify a transform to the reader so that the data will be stored in memory in a patient coordinate system. This insures the proper treatment of left/right and top/bottom in the resulting 3D models.
All the other parameters are self-explanatory except for the last. In this script, we know that the pipeline will only be executed once. To conserve memory, we set the GlobalReleaseDataFlag On. This allows the vtk pipeline execution to release data once it has been processed by a filter. For large medical datasets, this can mean the difference between being able to process a dataset and not.

set origin [expr ( $RESOLUTION / 2.0 ) * $PIXEL_SIZE * -1.0]
set minx [lindex $VOI 0]
set maxx [lindex $VOI 1]
set miny [lindex $VOI 2]
set maxy [lindex $VOI 3]
set minz [lindex $VOI 4]
set maxz [lindex $VOI 5]
#
# adjust y bounds for lower left coordinate system
#
set tmp $miny
set miny [expr $RESOLUTION - $maxy -1]
set maxy [expr $RESOLUTION - $tmp -1]
vtkImageReader reader;
    eval reader SetDataExtent 0 [expr $RESOLUTION - 1] 0 [expr $RESOLUTION - 1] $START_SLICE $END_SLICE 0 0
    reader SetDataVOI $minx $maxx $miny $maxy $minz $maxz
    reader SetFilePrefix $STUDY
    reader SetDataSpacing $PIXEL_SIZE $PIXEL_SIZE $SPACING
    eval reader SetDataOrigin $origin $origin [expr ( $START_SLICE - 1) * $SPACING]
    reader SetTransform $SLICE_ORDER
    reader SetHeaderSize $HEADER_SIZE
    reader SetDataMask 0x7fff;
    reader SetDataByteOrderTo$BYTE_ORDER
    reader SetDataScalarTypeToUnsignedShort
    [reader GetOutput] GlobalReleaseDataFlagOn;

Generate Triangles

The filter runs faster if we turn off gradient and normal calculations. Marching Cubes normally calculates vertex normals from the gradient of the volume data. In our pipeline, we will subsequently decimate the triangle mesh and smooth the resulting vertices. This processing invalidates the normals that can be calculated by Marching Cubes.

vtkMarchingCubes mcubes;
    mcubes SetInput [reader GetOutput]
    mcubes ComputeScalarsOff
    mcubes ComputeGradientsOff
    mcubes ComputeNormalsOff
    eval mcubes SetValue 0 $VALUE

Decimate The Triangles

There are often many more triangles generated by Marching Cubes than we need for rendering. Here we reduce the triangle count by eliminating triangle vertices that lie withinh a user-specified distance to the plane formed by neighboring vertices. We preserve an edges of triangles that are considered "features".

vtkDecimate decimator
    decimator SetInput [mcubes GetOutput]
    eval decimator SetInitialFeatureAngle $DECIMATE_ANGLE
    eval decimator SetMaximumIterations $DECIMATE_ITERATIONS
    decimator SetMaximumSubIterations 0
    decimator PreserveEdgesOn
    decimator SetMaximumError 1
    decimator SetTargetReduction $DECIMATE_REDUCTION
    eval decimator SetInitialError $DECIMATE_ERROR
    eval decimator SetErrorIncrement $DECIMATE_ERROR_INCREMENT

Smooth the Triangle Vertices

This filter uses LaPlacian smoothing to adjust triangle vertices as an "average" of neighboring vertices. Typically, the movement will be less that a voxel.

	
vtkSmoothPolyDataFilter smoother
    smoother SetInput [decimator GetOutput]
    eval smoother SetNumberOfIterations $SMOOTH_ITERATIONS
    eval smoother SetRelaxationFactor $SMOOTH_FACTOR
    eval smoother SetFeatureAngle $SMOOTH_ANGLE
    smoother FeatureEdgeSmoothingOff
    smoother BoundarySmoothingOff;
    smoother SetConvergence 0

Generate Normals

To generate smooth shaded models during rendering, we need normals at each vertex.

vtkPolyDataNormals normals
    normals SetInput [transformer GetOutput]
    eval normals SetFeatureAngle $FEATURE_ANGLE

Generate Triangle Strips

Triangle strips are a compact representation of large numbers of triangles. This filter processes our independent triangles before we write them to a file.

vtkStripper stripper
    stripper SetInput [smoother GetOutput]

Write the Triangles to a File

Finally, the last component of the pipeline write the triangles to a file.

vtkPolyDataWriter writer
    writer SetInput [stripper GetOutput]
    eval writer SetFilename $NAME.vtk
    writer SetFileType 2

Execute The Pipeline

All of the above script sets up connections in the pipeline. Up until this point, no files have been read nor triangles processed. This next section causes the pipeline to be executed. Normally, execution begins when the end of the pipeline, requests its input. Here, that would be the vtkPolyDataWriter. But for this script, we want to time the individual components of the pipeline, so we explicitly force a pipeline component execution with an Upadte method.

Set The Start Methods

These first few lines of tcl set the start methods for each component so we can see each filter start.

proc readerStart {} {global NAME; puts -nonewline "$NAME read took:\t"; flush stdout};
reader SetStartMethod readerStart
proc mcubesStart {} {global NAME; puts -nonewline "$NAME mcubes generated\t"; flush stdout};
proc mcubesEnd {} {
    global NAME
    puts -nonewline "[[mcubes GetOutput] GetNumberOfPolys]"
    puts -nonewline " polygons in "
    flush stdout
};
mcubes SetStartMethod mcubesStart
mcubes SetEndMethod mcubesEnd
proc decimatorStart {} {global NAME; puts -nonewline "$NAME decimator generated\t"; flush stdout};
proc decimatorEnd {} {
    global NAME
    puts -nonewline "[[decimator GetOutput] GetNumberOfPolys]"
    puts -nonewline " polygons in "
    flush stdout
};
decimator SetStartMethod decimatorStart
decimator SetEndMethod decimatorEnd
proc smootherStart {} {global NAME; puts -nonewline "$NAME smoother took:\t"; flush stdout};
smoother SetStartMethod smootherStart
proc writerStart {} {global NAME; puts -nonewline "$NAME writer took:\t"; flush stdout};
transformer SetStartMethod transformerStart
proc transformerStart {} {global NAME; puts -nonewline "$NAME transformer took:\t"; flush stdout};
stripper SetStartMethod stripperStart
proc stripperStart {} {global NAME; puts -nonewline "$NAME stripper took:\t"; flush stdout};

Execute Each Pipeline Filter

These next statements explicitly start each pipeline component. We use the tcl time command to show the elapsed time for each step in the process.

puts "[expr [lindex [time {reader Update;} 1] 0] / 1000000.0] seconds"
puts "[expr [lindex [time {mcubes Update;} 1] 0] / 1000000.0] seconds"
puts "[expr [lindex [time {decimator Update;} 1] 0] / 1000000.0] seconds"
puts "[expr [lindex [time {smoother Update;} 1] 0] / 1000000.0] seconds"
puts "[expr [lindex [time {normals Update;} 1] 0] / 1000000.0] seconds"
puts "[expr [lindex [time {stripper Update;} 1] 0] / 1000000.0] seconds"
writerStart
puts "[expr [lindex [time {writer Update;} 1] 0] / 1000000.0] seconds"
exit

And that's all there is to it. The models generated by the script can be viewed with renactors.tcl.

Questions / Comments

Bill Lorensen (lorensen@crd.ge.com)

Last Updated: Friday, 01-May-98 07:40:43