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

Extract A Volume Of Interest

Generate Triangles

Decimate The Triangles

Smooth The Triangle Vertices

Transform The Triangle Vertices

Generate Normals

Generate Triangles Strips

Write The Triangles

Execute The Pipeline

Set The Start Methods

Execute Each Pipeline Filter

References

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 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 SWAP] == 0} {set SWAP SwapBytesOff}
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 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 SMOOTH_EXPAND] == 0} {set SMOOTH_EXPAND -.34}
if {[info exists SMOOTH_SHRINK] == 0} {set SMOOTH_SHRINK .33}
if {[info exists VOI] == 0} {
    set VOI "0 [expr $RESOLUTION - 1] 0 [expr $RESOLUTION - 1] 0 $ZMAX]"
}
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). 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]
vtkVolume16Reader reader;
    eval reader SetDataDimensions $RESOLUTION $RESOLUTION
    eval reader SetFilePrefix $STUDY
    eval reader SetDataAspectRatio $PIXEL_SIZE $PIXEL_SIZE $SPACING
    eval reader SetDataOrigin $origin $origin [expr $START_SLICE * $SPACING]
    eval reader SetImageRange $START_SLICE $END_SLICE
    reader SetHeaderSize $HEADER_SIZE
    reader SetDataMask 0x7fff;
    eval reader $SWAP
    [reader GetOutput] GlobalReleaseDataFlagOn;

Extract a Volume of Interest

This is a new filter for vtk 1.2, allowing the selection of a volume of interest(VOI) and sampling rate for the volume. A VOI lets us select areas of interest, sometimes eliminating extraneous structures like the CT table bed. Subsampling often reduces an unmanagable dataset into a managable one. Be forewarned that subsampling can also eliminate important structures. It is best used to quickly select other pipeline parameters such as isosurface value.

vtkExtractVOI extractor
    extractor SetInput [reader GetOutput]
    eval extractor SetVOI $VOI
    eval extractor SetSampleRate $SAMPLE_RATE

Generate Triangles

This filter uses a number of features new to vtk 1.2. 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 [extractor 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 is new in vtk 1.1. It uses LaPlacian smoothing to adjust triangle vertices as an "average" of neighboring vertices. Typically, the movement will be less that a voxel. See Taubin [1] for a technical description of polyhedral smoothing. Check out Experiments with Smoothing for some more examples of this technique.

	
vtkSmoothPolyFilter smoother
    smoother SetInput [decimator GetOutput]
    eval smoother SetNumberOfIterations $SMOOTH_ITERATIONS
    eval smoother SetExpansionFactor $SMOOTH_EXPAND
    eval smoother SetContractionFactor $SMOOTH_SHRINK 
    eval smoother SetFeatureAngle $SMOOTH_ANGLE
    smoother FeatureEdgeSmoothingOff
    smoother BoundarySmoothingOff;
    smoother SetConvergence 0

Transform the Triangle Vertices

Here's a step that is often neglected. Medical images can be acquired in a variety of orders. For example, in CT, the data can be gathered from top to bottom (superior to inferior), or bottom to top (inferior to superior). In addition, MRI data can be acquired from left to right, right to left, front to back (anterior to posterior) or back to front. This filter will transform triangle vertices such that the resulting models will all "face" the viewer with a view up of (0,-1,0), looking down the positive z axis. Also, proper left-right correspondence will be maintained. That means the patient's left will always be left oin the generated models. Look in permutes.tcl defines the permutations and rotations for each order.

vtkTransformPolyFilter transformer
    eval transformer SetTransform $SLICE_ORDER
    transformer SetInput [smoother GetOutput]

Generate Normals

To generate smooth shaded models during rendering, we need normals at each vertex. This filter generates those normals and also optionally flips the normals. Normal flipping is required if the transformation caused a reversal of inside and outside. As in decimation, "features" can be retained.

vtkPolyNormals normals
    normals SetInput [transformer GetOutput]
    eval normals SetFeatureAngle $FEATURE_ANGLE
    if { $SLICE_ORDER == "si" || $SLICE_ORDER == "pa" || $SLICE_ORDER == "rl" } {
        normals FlipNormalsOn
    }

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 [normals GetOutput]

Write the Triangles to a File

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

vtkPolyWriter 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 vtkPolyWriter. 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 extractorStart {} {global NAME; puts -nonewline "$NAME extractor took:\t"; flush stdout};
extractor SetStartMethod extractorStart
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 normalsStart {} {global NAME; puts -nonewline "$NAME normals took:\t"; flush stdout};
normals SetStartMethod normalsStart
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 {extractor 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 {transformer 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.

References

Taubin, G., "A Signal Processing Approach to Fair Surface Design," Computer Graphics Proceedings, pp. 351-358, August 1995.

Questions / Comments

Bill Lorensen (lorensen@crd.ge.com)

Last Updated: Thursday, 30-Aug-01 17:41:07

GE Home Page | GE Research and Development