CS6630 - Scientific Visualization - Project 1

Carlos Eduardo Scheidegger
cscheid@cs.utah.edu
cs6630

Visualizing R^2->R^1 Data

Part 1: Height Fields

#!/usr/bin/env python
# CS 6630, Fall 2005, Project 1
# Carlos Eduardo Scheidegger
# convert dataset to vtk format

from os import sys

pts = map(lambda x: x.split(), file(sys.argv[1], 'r').readlines())
data = file(sys.argv[2], 'r').readlines()

print "# vtk DataFile Version 3.0"
print "vtk output"
print "ASCII"
print "DATASET POLYDATA"
print "POINTS %s float" % len(pts)

for (pt, datapt) in zip(pts, data):
    (ptx, pty, _) = pt
    print "%s %s %s" % (float(ptx), float(pty), float(datapt))

print "POINT_DATA %s" % len(pts)
convert.py: converts file to vtk
#!/usr/local/pkgs/vtk-4.2/bin/vtk
# Triangulate input, given as unstructured points
source "common.tcl"
package require vtk

proc triangulate { input output } {
    vtkDataSetReader reader
    reader SetFileName $input
    vtkDelaunay2D del
    del SetTolerance 0.01
    del SetAlpha 9
    del SetInput [reader GetOutput]
    vtkDataSetWriter writer
    writer SetFileName $output
    writer SetInput [del GetOutput]
    writer Update
}

about
set inputFileName  "$VTK_DATA/[lindex $argv 0]"
set outputFileName "$VTK_DATA/[lindex $argv 1]"
triangulate $inputFileName $outputFileName
exit
triangulate.tcl: computes and saves alpha shape
#!/usr/local/pkgs/vtk-4.2/bin/vtk
package require vtk
source "common.tcl"
about

[vtkPNMReader reader] SetFileName $VTK_DATA/MtHood.pgm
[vtkImageDataGeometryFilter geometry] SetInput [reader GetOutput]

vtkWarpScalar warp
warp SetInput [geometry GetOutput]
warp SetScaleFactor 0.28575 # explained in the report

vtkMergeFilter merge
merge SetGeometry [warp GetOutput]
merge SetScalars [reader GetOutput]

[vtkDataSetMapper mapper] SetInput [merge GetOutput]
[vtkActor actor] SetMapper mapper

[vtkRenderer renderer] AddActor actor
setColorByName renderer {SetBackground} $cobalt

vtkRenderWindow renderWindow
renderWindow AddRenderer renderer
renderWindow SetSize 512 512

vtkRenderWindowInteractor interactor
interactor SetRenderWindow renderWindow
interactor Initialize
interactor SetInteractorStyle [vtkInteractorStyleTrackballCamera _]
view.tcl: visualizes alpha shapes
#!/usr/local/pkgs/vtk-4.2/bin/vtk

package require vtk

source "common.tcl"

about

[vtkPNMReader reader] SetFileName $VTK_DATA/MtHood.pgm
[vtkImageDataGeometryFilter geometry] SetInput [reader GetOutput]

vtkWarpScalar warp
warp SetInput [geometry GetOutput]
warp SetScaleFactor 0.28575 # explained in the report

vtkMergeFilter merge
merge SetGeometry [warp GetOutput]
merge SetScalars [reader GetOutput]

[vtkDataSetMapper mapper] SetInput [merge GetOutput]
[vtkActor actor] SetMapper mapper

[vtkRenderer renderer] AddActor actor
setColorByName renderer {SetBackground} $cobalt

vtkRenderWindow renderWindow
renderWindow AddRenderer renderer
renderWindow SetSize 512 512

vtkRenderWindowInteractor interactor
interactor SetRenderWindow renderWindow
interactor Initialize
interactor SetInteractorStyle [vtkInteractorStyleTrackballCamera _]
count.tcl: visualizes Mt. Hood

In the first part of the project, we want to generate visualizations of scalar fields in R^2, by a technique known as height fields. With height fields, the scalar value of each point in R^2 is deformed in the z-direction. This way, the users can notice changes in the data by analyzing changes in the shape of the height field. Naturally, topographical maps are the most common example of data visualized in this fashion. We will be visualizing real topographical data, and simulation data, both provided by the instructing staff.

The first dataset to visualize was provided by the instructing staff, and we were asked to convert it to a format VTK understands. I did this in two steps: first, I converted the raw ASCII points to an ASCII VTK unstructured point set. I created a small script in python convert.py, shown on the right, to do this. The script acts as a filter, taking the raw ASCII data on stdin and outputting the VTK dataset on stdout.

Since this data consists only of a set of points, and we want to generate visualizations a height-field, we need to generate a triangulation of these. I used vtk to generate, specifically, one particular alpha shape of the point set [Edels83]. Alpha shapes are a subset of the Delaunay complex, and can be computed using the vtkDelaunay2D filter. The file triangulated.vtk is the alpha shape of the supplied point set with alpha=20. The vtk script triangulate.tcl used to generate and save the triangulation is also shown on the right. Delaunay triangulations are a good choice for generating triangles from a fixed set of points because they tend to generate triangles of good quality if the input points are well-spaced. Particularly, they maximize the minimum angle of all triangles from the triangulation [deBerg00].

The next step is to actually generate a visualization of the first dataset. For that, I created a simple vtk script that takes an unstructured point set, generates its alpha shape, and shows it on the screen. The script is named view.tcl and takes as input on the command line the filename and alpha value, and visualizes the resulting alpha shape.

This figure was generated by running view.tcl as follows:

[cscheid@trust]$ ./view.tcl triangulated.vtk 20
can't find package vtkHybridTCL 4.2
can't find package vtkParallelTCL 4.2
can't find package vtkPatentedTCL 4.2
CS6630 - Scientific Visualization
Prof. Chuck Hansen
---
Carlos Eduardo Scheidegger - cscheid@cs.utah.edu
---
We're running VTK version 4.2.1

Figure 1: visualizing the torso

The next heightmap visualization done was from a scalar field of the height of Mt. Hood. For this, I adapted the image warping example given by the instructors. The biggest change I've made to the script (aside from the scale, which will be explained later) was to replace vtkStructuredPointsGeometryFilter, which is deprecated in VTK 4.2, by vtkImageDataGeometryFilter. Also, since we're working with a PGM image, which is always grayscale, the vtkImageLuminance filter is not necessary. I then simply used the scalars as they were given.

To get the right scaling factor (and then, visualize the correct height for Mt. Hood), I first found out Mt. Hood's elevation on Wikipedia (3.429km). Then, I assumed that the base of the image spanned all of the 3.429km, since Mt. Hood is in Portland, OR, and Portland is at sea level (this is obviously not strictly true, but seems good enough). Then, I looked at a map of Mt. Hood at Google Maps and eyeballed an estimate of the size of Mt. Hood at around 12km. This gives a necessary warp of around 0.28 so that the height appears the same. The result is shown in the following figure. 

[cscheid@trust]$ ./mount.tcl
can't find package vtkHybridTCL 4.2
can't find package vtkParallelTCL 4.2
can't find package vtkPatentedTCL 4.2
CS6630 - Scientific Visualization
Prof. Chuck Hansen
---
Carlos Eduardo Scheidegger - cscheid@cs.utah.edu
---
We're running VTK version 4.2.1

Figure 2: visualizing Mt. Hood

Part 2: Contour Maps

Another strategy for visualizing scalar fields (and one that is more easily extensible to higher dimensions) is the use of contour maps. Contour maps display only a part of the data at a time. Particularly, they display the preimage of a certain scalar in the dataset at each time. In 2D, the preimages of a continuous scalar fields are lines, and in 3D, they are surfaces. Contour maps can be very effective for visualizing quantitative data, because the user is allowed to pick a certain measurement value and investigate the shape of the manifold that has this value. Additionally, in closed compact spaces without boundary, the manifolds of the contour maps are themselves without boundary. In 2D, this means that every line forms a loop. In 3D, it means that there are on boundaries on the resulting 2D manifold. Since most datasets (not all - datasets with toroidal topology can be without boundary) are compact but with boundary, this may not be the case. Then, the contour lines that go to the boundary might have boundaries themselves.

To experiment with contour maps, I implemented a vtk/tk application in which several 2D datasets can be investigated. The two provided datasets are an infrared image of the human body and a (I assume) MRI scan of a human brain. I provided some datasets of my own, describing distance fields to lines in 2D. 

if { [catch {set VTK_TCL $env(VTK_TCL)}] != 0} { set VTK_TCL ".." }
if { [catch {set VTK_DATA $env(VTK_DATA)}] != 0} { set VTK_DATA "../vtkdata" }

proc invalidate { } {
    global currentActive
    global renWin

    set mnx 0
    set mxx 0
    set mny 0
    set mxy 0
    set mnz 0
    set mxz 0

    scan [[reader$currentActive GetOutput] GetBounds] "%f %f %f %f %f %f" mnx mxx mny mxy mnz mxz
    set mpx [expr [expr $mxx + $mnx] / 2.0]
    set mpy [expr [expr $mxy + $mny] / 2.0]
    set mpz [expr [expr $mxz + $mnz] / 2.0]

    set camera [ren GetActiveCamera]
    $camera ParallelProjectionOn
    $camera SetPosition $mpx $mpy [expr $mnz - 1]
    $camera SetFocalPoint $mpx $mpy 1
    $camera SetParallelScale [expr [expr $mxy - $mny] / 2.0]

    $renWin Render
}

proc changeImage { v } {
    global currentActive
    ren RemoveActor actor$currentActive
    set currentActive $v
    ren AddActor actor$currentActive
    invalidate
}

proc changeIsoValue { im v } {
    contour$im SetValue 0 $v
    invalidate
}

proc loadDataSet { fileName name value caption radiobuttonframe sliderframe rbvar rbcommand slcommand } {
    global VTK_DATA
    vtkDataSetReader reader$name
    reader$name SetFileName "$VTK_DATA/$fileName"
    reader$name Update
    set mins 0
    set maxs 0
    scan [[reader$name GetOutput] GetScalarRange] "%f %f" mins maxs
    set res [expr [expr $maxs - $mins] / 100]
    radiobutton $radiobuttonframe.$name -text $caption \
           -value $value -command $rbcommand -variable $rbvar
    pack $radiobuttonframe.$name -side top -padx 2 \
           -pady 2 -expand 0 -fill none -anchor w
    frame $sliderframe.$name
    label $sliderframe.$name.label -text $caption
    pack $sliderframe.$name.label -side left -padx 2  \
           -pady 2 -expand 0 -fill none -anchor w
    scale $sliderframe.$name.slider -from $mins -to $maxs \
	-resolution $res -orient horizontal -command $slcommand
    pack $sliderframe.$name.slider -side right -padx 2 \
          -pady 2 -expand 0 -fill none -anchor w
    pack $sliderframe.$name -side top -padx 0 -pady 2 \
         -expand 0 -fill none -anchor w
    [vtkContourFilter contour$name] SetValue 0 $mins
    contour$name SetInputConnection [reader$name GetOutputPort]
    [vtkPolyDataMapper polydatamapper$name] SetInput [contour$name GetOutput]
    [vtkActor actor$name] SetMapper polydatamapper$name
    contour$name Update
}

set datasetsLoadedRegister 0

proc smartLoadDataSet { fileName caption } {
    global datasetsLoadedRegister
    set changeRBProc [list changeImage $datasetsLoadedRegister ]
    set changeSliderProc [list changeIsoValue $datasetsLoadedRegister ]

    loadDataSet $fileName $datasetsLoadedRegister \
	$datasetsLoadedRegister $caption \
	{.top.left.radiobuttons} {.top.left.sliders} rbValue \
	$changeRBProc $changeSliderProc

    set datasetsLoadedRegister [expr $datasetsLoadedRegister + 1]
}


source "common.tcl"

### Create the top level stuff and basic frames

toplevel .top -visual best
wm title .top "Contouring 2D Data"

frame .top.left -width 256 -borderwidth 2 -relief ridge
frame .top.right -width 512 -height 512 -borderwidth 2 -relief ridge

pack .top.left .top.right -side left -expand 1 -fill both

### Create the render window

vtkTkRenderWidget .top.right.ren -width 512 -height 512
BindTkRenderWidget .top.right.ren

set renWin [.top.right.ren GetRenderWindow]
vtkRenderer ren
$renWin AddRenderer ren

# Test the render window
label .top.left.label -text "Control Options"
pack .top.left.label -side top -expand 0

set currentActive 0
set rbValue 0

# Set radiobutton and slider frames
frame .top.left.radiobuttons -borderwidth 2 -relief ridge
frame .top.left.sliders -borderwidth 2 -relief ridge

smartLoadDataSet "brain.vtk" "Brain"
smartLoadDataSet "body.vtk" "Body"

for {set i 0} {$i<$argc} {set i [expr $i+2]} {
    smartLoadDataSet [lindex $argv $i] [lindex $argv [expr $i+1]]
}

pack .top.left.radiobuttons -expand 0 -side top -pady 20
pack .top.left.sliders -expand 0 -side top -pady 20

# Final setup
changeImage 0
vtkRenderWindowInteractor iren
iren SetRenderWindow $renWin
iren SetInteractorStyle [vtkInteractorStyleTrackballCamera _]
contour.tcl: VTK/TK contouring app

My application uses the VolTkInteractor package provided by the instructing staff. Technically, the TkInteractor would be more appropriate, but I couldn't get it to work (it requires a vtkScalars vtk class that neither VTK 4.2 nor 4.5 seem to have). This package allows me to integrate the VTK window inside the tk widget tree. The entire application then consists of a single window, as follows:


Figure 3: Main window for contouring application
The correct ranges of the sliders are determined on a per-dataset basis, by using the GetScalarRange method that is provided by vtkDataSetSource objects. The correct camera positioning is determined likewise, this time by using the GetBounds method. Instead of simply using ResetCamera (which resets the camera to display the entire dataset and will automatically zoom to small contours, which might confuse the users), I determine the maximum bounds with GetBounds and then use SetPosition, SetFocalPoint and SetParallelScale manually. My tk script allows for an arbitrary number of datasets. The two default datasets used are always included, and the remaining ones are given from the command-line, by giving the vtk dataset name and a caption for the radio button and slider. Internally, they are handled very uglily (or cleverly, take your pick) by generating different variable names in runtime. Each dataset is given a unique identifier, and widgets, datasets, etc. are created by concatenating this unique dataset to fixed variable names. This should have been much cleaner, by using arrays or maybe classes from [incr tcl], but I haven't done so. In the following, there is a figure of my script running with five different datasets:

Figure 4: Many datasets
Although contouring is a very effective visualization tool, there are some requirements on the dataset for contouring to make sense. First of all, the scalar field should be continuous. When this is the case, contour lines are always nested (this can be proven by Stokes' theorem). This allows navigating the data sensibly, since the user can predict locally what will happen when pushing the slider either way. The other desirable property is monotonicity. By this, I mean that increasing scalar values should indicate an increase of some sensible quantity in the real world. If this "conceptual" monotonicity is broken, then the user will not be able to tell whether he/she should push the slider up or down to get the desired effect. Additionally, if the scalar field does not have this monotonicity, there might not even be a bijection between the scalar values and the physical quantities, making effective visualization very hard. Most observational scalar fields will be created so that they respect these properties. For example, most topographical datasets respect this (an example of a place where it wouldn't is any cliff where there are "two floors": two points with the same latitude and longitude but with different heights. Whichever point one picks for the topo map will incur in a discontinuity somewhere). Most scalar fields that directly represent some physical quantity will also repsect these (examples where this might not be the case are shock data, for example from explosions or supersonic flows) 
#!/usr/bin/env python

from math import *
from random import *

f = file('/dev/stdin', 'r')

lines = f.readlines()

(width, height) = lines[2].split()

width = int(width)
height = int(height)
counter = 0

points = []

def dist(pt1, pt2):
    (x1,y1) = pt1
    (x2,y2) = pt2
    dx = x2-x1
    dy = y2-y1
    return sqrt(dx*dx+dy*dy)

def dist_to_shape(pt):
    mn = dist(pt, points[0])
    for pt2 in points[1:]:
        d = dist(pt, pt2)
        if d < mn:
            mn = d
    return mn

for x in xrange(width):
    for y in xrange(height):
        if (int(lines[(y * width + x) * 3+4]) == 0):
            points.append( (float(2*x)/float(width)-1.0,
                            float(2*y)/float(height)-1.0) )

print "# vtk DataFile Version 2.0"
print "vtk output"
print "ASCII"
print "DATASET STRUCTURED_POINTS"
print "DIMENSIONS %s %s 1" % (width, height)
print "ORIGIN 0 0 0"
print "SPACING 1 1 1"
print "POINT_DATA %s" % (width * height)
print "SCALARS values float"
print "LOOKUP_TABLE default"

for x in xrange(width):
    for y in xrange(height):
        ptx = uniform(-1, 1)
        pty = uniform(-1, 1)
        print "%s" % int(256 * dist_to_shape( (float(2*x)/float(width)-1.0,
                                               float(2*y)/float(height)-1.0) ))
distance.py: computing distance fields
My additional datasets created for the contouring are distance fields. I find distance fields interesting because contours of distance fields create what are called offset surfaces, which are geometrical entities with many applications, including computer-aided design and motion planning. The interesting thing about offset surfaces that I wanted to visualize is the fact that simply-connected surfaces might generated non-simply connected offset surfaces. This is important for motion planning because it might create "isolated islands" between which movement is not possible. To create the distance fields, I coded a simple python script that read a .ppm image, in which a black point denoted a point on the surface. Then I stored all the surface points in a list, and for every point in the image, I computed the distance to the surface by simply iterating through all points. This is very very inefficient (worst-case O(n^4)), but since I didn't care about preprocessing time, it was good enough for my purposes. I'm including in my submission two datasets that I created, both featuring the phenomenon I described above. I'm also including the python script. To run my application, you can either run it directly from the command line, or use the driver scripts contour_driver.sh and contour_driver_2.sh.

Figure 5: Offset surfaces as contours of distance fields

Figure 6: A different offset surface

References

[1] H. Edelsbrunner, D.G. Kirkpatrick and R. Seidel. On the shape of a set of points in the plane. IEEE Trans. Inform. Theory 29, 551-559 (1983).

[2] M. de Berg et al. Computational Geometry - Algorithms and Applications. Springer-Verlag, 2000.