Using pick to create curves

VisIt provides a 2D and 3D lineout capabilities that allow you to extract data along a line. What about a more complex curve though? It is possible to extract data along a curve using VisIt's Pick capability and some Python scripting. The following code example shows how to create a curve based on samples extracted from the surface of a cylinder. A path is created along a circle in 3D and the script picks at each sample point in the path and produces a curve file from that data to produce the Lineout curve.

The script

import math, random

#
# Vector functions
#
def vec_normalize(vec):
    R = math.sqrt(vec[0]*vec[0] + vec[1]*vec[1] + vec[2]*vec[2])
    return (vec[0] / R, vec[1] / R, vec[2] / R)

def vec_add(a, b):
    return (a[0] + b[0], a[1] + b[1], a[2] + b[2])

def vec_mult_float(vec, f):
    return (vec[0] * f, vec[1] * f, vec[2] * f)

def vec_cross(a, b):
    return (a[1]*b[2] - a[2]*b[1],
            a[2]*b[0] - a[0]*b[2],
            a[0]*b[1] - a[1]*b[0])

#
# Write the path to a file and plot it.
#
def DebugLocations(locations):
    f = open("points.3D", "wt")
    f.write("X\n");
    f.write("Y\n");
    f.write("Z\n");
    f.write("Index\n");
    for index,loc in enumerate(locations):
        f.write("%g %g %g %d\n" % (loc[0], loc[1], loc[2], index))
    f.close()
    OpenDatabase("points.3D")
    AddPlot("Pseudocolor", "Index")
    DrawPlots()
    pc = PseudocolorAttributes()
    pc.pointSizePixels = 20
    pc.pointType = pc.Sphere
    SetPlotOptions(pc)

#
# Splits a string into a list of strings.
#
def SplitLines(s):
    lines = s.split("\n")
    return [line for line in lines if line != '']

#
# Extracts values from VisIt's pick output.
#
def ExtractPickValue(s):
    lines = SplitLines(s)
    retval = 0.
    print lines
    if "<nodal>" in lines[1]:
        sum = 0.
        nSum = 0
        # Extract all of the values in the pick and average them. Not
        # exactly what we want but it's something.
        for line in lines[2:]:
            e = line.find("=")
            if e != -1:
                tok = line[e+2:]
                value = float(tok)
                sum += value
                nSum += 1
        retval = sum / float(nSum)
    elif "<zonal>" in lines[1]:
        e = lines[1].find("=")
        if e != -1:
            tok = lines[1][e+2:]
            retval = float(tok)
    return retval

#
# Calculates the sample points for a circle in 3D
#
# origin : A tuple of 3 floats containing the circle's origin
# normal : A tuple of 3 floats containing the normal of the plane
#          that contains the circle
# radius : The circle's radius
# npts   : The number of samples around the circle
#
def CircularPickLocationsAboutAxis(origin, N, radius, npts = 10):
   normal = vec_normalize(N)
   a2 = (0., 1., 0.)
   up = vec_normalize(vec_cross(a2, normal))
   right = vec_normalize(vec_cross(up, normal))
   locations = []
   for i in range(npts):
       angle = math.pi * 2. * float(i) / float(npts-1)
       xcoeff = radius * math.cos(angle)
       ycoeff = radius * math.sin(angle)
       loc = vec_add(origin,vec_add(vec_mult_float(right, xcoeff),vec_mult_float(up,ycoeff)))
       locations.append(loc)
   return locations


#
# Extracts data along a path and returns a tuple
#
# path : A tuple of 3-float tuples containing the 3D sample points
#
def ExtractDataAlongPath(path):
    p = GetPickAttributes()
    p.variables = ("default")
    p.displayIncidentElements = 0
    p.showNodeId = 0
    p.showNodeDomainLogicalCoords = 0
    p.showNodeBlockLogicalCoords = 0
    p.showNodePhysicalCoords = 0
    p.showZoneId = 0
    p.showZoneDomainLogicalCoords = 0
    p.showZoneBlockLogicalCoords = 0
    p.doTimeCurve = 0
    p.conciseOutput = 1
    p.showTimeStep = 0
    p.showMeshName = 0
    p.blockPieceName = "domain"
    p.groupPieceName = "block"
    p.displayGlobalIds = 0
    p.displayPickLetter = 0  # Change this to 1 if you want to see the pick points
    p.meshCoordType = 0
    p.createSpreadsheet = 0
    p.floatFormat = "%g"
    p.timePreserveCoord = 0
    SetPickAttributes(p)
   
    # Write the points to a file for debugging.
    #DebugLocations(path)

    SuppressQueryOutputOn()
    samples = []
    for pt in path:
        ZonePick(pt, ("default"))
        samples.append(ExtractPickValue(GetPickOutput()))
    SuppressQueryOutputOff()
    return samples

#
# Calculates a lineout along a path using pick and plots the 
# lineout curve in a new window.
#
# path : A tuple of 3-float tuples containing the 3D sample points
#
def LineoutAlongPath(path):
    samples = ExtractDataAlongPath(path)

    # Write the samples to a curve file.
    rindex = random.randrange(1000)
    dbName = "samples%04d.curve" % rindex
    f = open(dbName, "wt")
    f.write("# samples\n")
    for index,s in enumerate(samples):
        f.write("%d %g\n" % (index, samples[index]))
    f.close()

    # Plot the curve
    AddWindow()
    DeleteAllPlots()
    OpenDatabase(dbName)
    AddPlot("Curve", "samples")
    DrawPlots()

#
# Performs a circular lineout in 3D
#
# origin : A tuple of 3 floats containing the circle's origin
# normal : A tuple of 3 floats containing the normal of the plane
#          that contains the circle
# radius : The circle's radius
# npts   : The number of samples around the circle
#
def CircularLineout(origin, normal, radius, npts = 10):
    path = CircularPickLocationsAboutAxis(origin, normal, radius, npts)
    LineoutAlongPath(path)

#
# Demonstrates the circular lineout.
#
def demo():
    OpenDatabase("/usr/gapps/visit/data/noise.silo")
    AddPlot("Pseudocolor", "airVf")
    DrawPlots()

    p1 = (-6, -6, -6)
    p2 = (6, 6, 6)
    AddOperator("Cylinder")
    cyl = CylinderAttributes()
    cyl.point1 = p1
    cyl.point2 = p2
    cyl.radius = 4.
    SetOperatorOptions(cyl)

    # 83% of the distance from p1 to p2 is where our circle's origin will be
    t = 0.83

    # Do a circular lineout around the edges of the cylinder
    origin = vec_add(vec_mult_float(p1, 1.-t), vec_mult_float(p2, t))
    normal = vec_normalize(vec_add(p2, vec_mult_float(p1, -1.)))
    CircularLineout(origin, normal, cyl.radius, npts = 50)

The results

Picking on cylinder surface normals

Picking on cylinder surface normals as created by the DeferExpression operator is not as straightforward as picking on scalar values. For one thing, pick does not seem capable of returning data values because the expression value calculation has been deferred until after pick. Another problem is that picking on vector plots is not terribly easy since only picks on the glyphs return values.

Here is a possible workaround script

• Create the surface normals and export vector components to a new file
• Pick on the vector components separately to avoid the problems inherent in the pick on Vector plots

# All of the above script is needed to run this script...

def pick_surface_normals(db, p1, p2, rad, t):
    import os
    OpenDatabase(db)

    # Define the surface normal expression that we'll use to plot vectors
    DefineVectorExpression("sn", "-point_surface_normal(Mesh)")
    DefineScalarExpression("sn0", "sn[0]")
    DefineScalarExpression("sn1", "sn[1]")
    DefineScalarExpression("sn2", "sn[2]")

    AddPlot("Pseudocolor", "sn0")
    AddOperator("Cylinder")
    cyl = CylinderAttributes()
    cyl.point1 = p1
    cyl.point2 = p2
    cyl.radius = rad
    SetOperatorOptions(cyl)

    # You must have DeferExpression activated in your settings since it
    # is not enabled by default
    AddOperator("DeferExpression")
    de = DeferExpressionAttributes()
    de.exprs =("sn0", "sn1", "sn2", "sn")
    SetOperatorOptions(de)

    DrawPlots()

    vectorComps = ("sn0", "sn1", "sn2")

    # Export the database - it's the only way to get past pick
    # on DeferExpression. This could result in many VTK files if
    # you're exporting a domain-decomposed database.
    dbAtts = ExportDBAttributes()
    dbAtts.db_type = "VTK"
    dbAtts.filename = "sn"
    dbAtts.variables = vectorComps
    opts = GetExportOptions("VTK")
    opts['Binary format'] = 1
    ExportDatabase(dbAtts, opts)

    # Clean up after the export since we want to plot the exported
    # database from now on.
    DeleteAllPlots()

    # Undefine the expressions since they'll conflict with the new db.
    DeleteExpression("sn0")
    DeleteExpression("sn1")
    DeleteExpression("sn2")

    # Figure out the path along the outside of the cylinder
    origin = vec_add(vec_mult_float(p1, 1.-t), vec_mult_float(p2, t))
    normal = vec_normalize(vec_add(p2, vec_mult_float(p1, -1.)))
    path = CircularPickLocationsAboutAxis(origin, normal, cyl.radius, npts = 50)

    # See which file was saved by the export.
    db2 = "sn.visit"
    try:
        s = os.stat(db2)
    except OSError:
        db2 = "sn.vtk"

    # Add a new plot
    OpenDatabase(db2)
    AddPlot("Pseudocolor", "sn0")
    DrawPlots()
    sn0 = ExtractDataAlongPath(path)
    ChangeActivePlotsVar("sn1")
    sn1 = ExtractDataAlongPath(path)
    ChangeActivePlotsVar("sn2")
    sn2 = ExtractDataAlongPath(path)
    DeleteAllPlots()

    # Combine to make a list of (sn0, sn1, sn2) tuples
    sn = zip(sn0, sn1, sn2)
    return (path, sn)

#
# Save the points and vec to an easy file format.
#
def SaveXmdv(filename, pts, vec):
    f = open(filename, "wt")
    f.write("6 %d 12\n" % len(pts))
    f.write("x\ny\nz\nsnx\nsny\nsnz\n")
    # Write some bogus extents
    for i in range(6):
        f.write("0.000000        1.000000        10\n")
    # Write columns of data values
    for pt,sn in zip(pts, vec):
        f.write("%1.7f %1.7f %1.7f %1.7f %1.7f %1.7f\n" % (pt[0],pt[1],pt[2],sn[0],sn[1],sn[2]))
    f.close()

def demo():
    p1 = (-6, -6, -6)
    p2 = (6, 6, 6)
    rad = 4.
    t = 0.83

    # Extract the data
    db = "/usr/gapps/visit/data/noise.silo"
    path_sn = pick_surface_normals(db, p1, p2, rad, t)

    # Let's do a check by plotting the surface normal vectors.
    path = path_sn[0]
    sn = path_sn[1]
    SaveXmdv("normals.okc", path, sn)

    # Do a plot of the original database with the Cylinder operator
    OpenDatabase(db)
    AddPlot("Pseudocolor", "airVf")
    AddOperator("Cylinder")
    cyl = CylinderAttributes()
    cyl.point1 = p1
    cyl.point2 = p2
    cyl.radius = rad
    SetOperatorOptions(cyl)
    DrawPlots()

    OpenDatabase("normals.okc")
    DefineVectorExpression("sn", "{snx, sny, snz}")
    AddPlot("Vector", "sn")
    DrawPlots()

    return

Surface normal picking results