Creating curves from contours

Overview

The lineout operation in VisIt plots plots values versus distance along a straight line. As of this writing (August 26, 2011), there is no convenient way to perform a lineout along a curve. However, contours lines created using the Isosurface operator (with the isosurface variable set to define the shape of the curve) and the Pseudocolor plot (with the plot variable set to the values you wish to extract) gets us very close, with values and line segments. All that remains is to "unwrap" or "flatten" each of the contours into a linear sequence. The python script below accomplishes that, by joining individual line segments into multi-segment "polylines" and writing their values to a .curve file as a function of accumulated distance along the polyline.

To use this routine:

  • Create a Pseudocolor plot of your value variable
  • Add an Isosurface operator and set it to use your level set variable
  • In File -> Set save options..., set your output file type to VTK and ensure "Binary" is NOT checked
  • Run this script; its first argument is the name of this file you just saved, and its second is an output .curve file name
  • Open the resulting .curve file in VisIt

Source code

# File:       vtklines_to_curves.py
# Programmer: Jeremy Meredith
# Date:       August 26, 2011
#
# Description:
#   Input a file containing a VTK_POLYDATA with only 2-point line segments
#   and either a point or cell scalar.  Coalesce the short line segments into
#   polylines.  Output each polyline as a function of scalar-value vs distance.
#
# Usage:
#   python vtklines_to_curves.py <inputfile.vtk> <outputfile.curve>
#
import sys, os, math


# ----------------------------------------------------------------------------
#  Parse the points, line segments and cell or point scalars from a VTK file
# ----------------------------------------------------------------------------
def ParseFile(fn):
    f = open(fn)
    try:
        npts = 0
        pts = []
        nlines = 0
        lines = []
        vals  = []
        name  = "values"
        ispoint = -1
        # modes: 0 = searching for POINTS
        #        1 = reading points
        #        2 = looking for LINES
        #        3 = reading lines
        #        4 = looking for CELL_DATA, POINT_DATA, or SCALAR
        #        5 = found scalar, skipping LOOKUP_TABLE
        #        6 = reading scalar values
        mode = 0
        for l in f:
            words = l.split()
            if len(words) == 0:
                None
            elif mode == 0:
                if words[0] == "POINTS":
                    npts = int(words[1])
                    if npts == 0:
                        raise "Couldn't extract number of points"
                    mode = 1
            elif mode == 1:
                for w in words:
                    pts.append(float(w))
                if len(pts) >= npts*3:
                    mode = 2
            elif mode == 2:
                if words[0] == "LINES":
                    nlines = int(words[1])
                    if nlines == 0:
                        raise "Couldn't extract number of lines"
                    nints  =  int(words[2])
                    if nints != 3*nlines:
                        raise "Expected only 2-point line segments"
                    mode = 3
            elif mode == 3:
                if int(words[0]) != 2:
                    raise "Expected only 2-point line segments"
                seg = [int(words[1]),int(words[2])]
                lines.append(seg)
                if len(lines) >= nlines:
                    mode = 4
            elif mode == 4:
                if words[0] == "CELL_DATA":
                    ispoint = 0
                elif words[0] == "POINT_DATA":
                    ispoint = 1
                elif words[0] == "SCALARS":
                    if ispoint < 0:
                        raise "Found scalars before CELL_DATA or POINT_DATA"
                    name = words[1]
                    mode = 5
            elif mode == 5:
                if words[0] == "LOOKUP_TABLE":
                    mode = 6
            elif mode == 6:
                for w in words:
                    vals.append(float(w))
                if ((ispoint==1 and len(vals) >= npts) or
                    (ispoint==0 and len(vals) >= nlines)):
                    mode = 7
            else: # mode == 7
                break
    finally:
        f.close()
    return (name, ispoint, pts, lines, vals)

# ----------------------------------------------------------------------------
#  Operations on points
# ----------------------------------------------------------------------------
def calcdist(pt0, pt1):
    dx = pt0[0] - pt1[0]
    dy = pt0[1] - pt1[1]
    dz = pt0[2] - pt1[2]
    dist = math.sqrt(dx*dx + dy*dy + dz*dz)
    return dist

def calcmid(pt0, pt1):
    x = (pt0[0] + pt1[0]) / 2.
    y = (pt0[1] + pt1[1]) / 2.
    z = (pt0[2] + pt1[2]) / 2.
    return (x,y,z)

def getpoint(pts, i):
    return (pts[3*i + 0], pts[3*i + 1], pts[3*i + 2])



# ----------------------------------------------------------------------------
#  Main routine
# ----------------------------------------------------------------------------

# get command-line arguments
if len(sys.argv) != 3:
    print "Usage:",sys.argv[0],"<inputfile.vtk> <outputfile.curve>"
    sys.exit(1)

infile = sys.argv[1]
if not os.path.exists(infile):
    print "Error:",infile,"didn't exist."
    sys.exit(1)

# parse the file
(name, ispoint, pts, lines, vals) = ParseFile(infile)

# for cell vars, put the point vals in a dictionary
linevals = {}
if not ispoint:
    for i in range(len(lines)):
        seg = lines[i]
        linevals[seg[0],seg[1]] = vals[i]
        linevals[seg[1],seg[0]] = vals[i]

# group the line segments into polylines, first pass
polylines = []
for seg in lines:
    found = False
    for i in range(len(polylines)):
        # if we can stick it on either end of a polyline, do it
        if seg[0] == polylines[i][0]:
            polylines[i].insert(0,seg[1])
            found = True
        elif seg[1] == polylines[i][0]:
            polylines[i].insert(0,seg[0])
            found = True
        elif seg[0] == polylines[i][-1]:
            polylines[i].append(seg[1])
            found = True
        elif seg[1] == polylines[i][-1]:
            polylines[i].append(seg[0])
            found = True
        if found:
            break
    if found == False:
        polylines.append(seg)

# now group up the polylines until we can't do it anymore
done = False
while not done:
    done = True
    for i in range(len(polylines)):
        for j in range(i):
            polyi = polylines[i]
            polyj = polylines[j]
            # if either end of the polylines match, merge them
            if polyi[0] == polyj[0]:
                polylines.remove(polyj)
                polyi.reverse()
                polyi.extend(polyj[1:])
                done = False
            elif polyi[-1] == polyj[0]:
                polylines.remove(polyj)
                polyi.extend(polyj[1:])
                done = False
            elif polyi[0] == polyj[-1]:
                polylines.remove(polyi)
                polyj.extend(polyi[1:])
                done = False
            elif polyi[-1] == polyj[-1]:
                polylines.remove(polyj)
                polyj.reverse()
                polyi.extend(polyj[1:])
                done = False

            if done == False:
                break
        if done == False:
            break

# now write the curve file
outfile = sys.argv[2]
out = open(outfile, 'w')
for p in range(len(polylines)):
    polyline = polylines[p]
    # note it's different between cell- and point-centered values
    out.write("# %s%04d\n" % (name, p))
    if ispoint:
        pos = 0.0
        for i in range(len(polyline)-1):
            if i == 0:
                out.write("0 %f\n" % vals[polyline[0]])
            pt0 = getpoint(pts, polyline[i])
            pt1 = getpoint(pts, polyline[i+1])
            pos += calcdist(pt0, pt1)
            out.write("%f %f\n" % (pos, vals[polyline[i]]))
    else:
        pos = 0.0
        lastpt = (0,0,0)
        for i in range(len(polyline)-1):
            pt0 = getpoint(pts, polyline[i])
            pt1 = getpoint(pts, polyline[i+1])
            pt = calcmid(pt0,pt1)
            val = linevals[polyline[i],polyline[i+1]]
            if i == 0:
                out.write("0 %f\n" % val)
                lastpt = pt
            else:
                pos += calcdist(lastpt, pt)
                lastpt = pt
                out.write("%f %f\n" % (pos, val))
out.close()

Example

  • Open noise.silo
  • Create a Pseudocolor plot of "grad_magnitude"
  • Add a Slice operator, along the X-axis, with an intercept at x=2
  • Add an Isosurface operator, with its variable set to "hardyglobal" and Value set to "3.4"
  • Note that the plot has two contours: one is short, closed with a high maximum near 0.7, and the other is large, open, with a lower maximum under 0.5.

Lineout curve setup.png

  • Save the window as data.vtk
  • Run the above script: python vtklines_to_curves.py data.vtk lineout.curve
  • Open lineout.curve
  • Plot grad_magnitude0000 and grad_magnitude0001
  • This shows the above contour values of grad_magnitude versus distance along the curves.
    • Note, as above, one is short with a high maximum value, and the other is long with a lower maximum value.

Lineout curve result.png

Notes

Ideally, a new operator or other intrinsic VisIt method will replace this method eventually.

Note that the script assumes you have an ASCII VTK polydata file with one 2-segment lines, and only a single scalar variable. This is a safe assumption when starting with a Pseudocolor plot of a Isosurface operator saved to an ASCII VTK file using Save Window in VisIt.

The above script could be easily modified to work within VisIt by saving the current window to a filename and passing that directly to the Parse routine instead of using a value from sys.argv, and using a temporary file to write the result to, subsequently opening it for the user.

Also, note that Using pick to create curves presents a good alternative, and would simply require that you get your isocontour locations in sorted order via some other method, although as it performs pick on these locations may be less efficient.