Rectout

VisIt has lineouts which sample values along a line and produces curves. This script produces a "rectout" which is like a lineout except that it has width. The path along the line is clipped to produce an area dA over which the plotted field is scaled according to its current/original area and then the value is integrated to produce a value. The integrated value is plotted vs the distance along the line and plotted as a curve.

################################################################################
# rectout.py -- This script defines a function called Rectout that attempts to
#               perform an integration of rectangular boxes along a line and
#               save the results to a curve file that then gets visualized.
#               For each dA in the clipped rectangular region, we weight the
#               plotted variable by how much of the original cell that is present
#               before integrating.
#
# Usage:
#     Rectout(p0, p1, width, nbins)
#
#     p0    : A 3-tuple that specifies the starting point for the lineout
#     p1    : A 3-tuple that specifies the end point for the lineout
#     width : How wide the rectangle should be. The rectangle extends width/2 
#             units on either side of the lineout line.
#     nbins : The number of dA to sample along the line.
#
# Example:
#     Rectout((0.1666, 0.8333, 0), (0.8, 0.2, 0.), 0.25, 20)
#
# Note: Before using the Rectout function, set up a plot of a scalar variable.
#       Also, you must enable the DeferExpression operator before using.
#
#       It might be possible to make a Macro function that applies the Rectout
#       function to any plots that have a Lineout operator to make getting the
#       p0, p1 points more automatic.
#
# Programmer: Brad Whitlock
# Date: Tue Sep  8 16:49:22 PDT 2009
#
# Modifications:
#
################################################################################
import math

#
# Vector functions
#
def vec_magnitude(vec):
    return math.sqrt(vec[0]*vec[0] + vec[1]*vec[1] + vec[2]*vec[2])

def vec_normalize(vec):
    R = vec_magnitude(vec)
    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])

def GetRectoutClipPlanes(p0, p1, width, startT=0, endT=1):
    p0p1 = (p1[0] - p0[0], p1[1] - p0[1], p1[2] - p0[2])
    dist = vec_magnitude(p0p1)
    p0p1_norm = vec_normalize(p0p1)

    v1 = vec_cross((0,0,1), p0p1_norm)
    v2 = vec_mult_float(v1, -1)

    s0 = vec_add(vec_mult_float(p0, 1.-startT), vec_mult_float(p1, startT))
    s1 = vec_add(vec_mult_float(p0, 1.-endT), vec_mult_float(p1, endT))

    # assemble origin/normal pairs
    plane0 = (vec_add(s0, vec_mult_float(v1, width/2)), v2)
    plane1 = (vec_add(s0, vec_mult_float(v2, width/2)), v1)

    plane2 = (s0, p0p1_norm)
    plane3 = (s1, vec_mult_float(p0p1_norm, -1.))

    return (plane0, plane1, plane2, plane3)


def RectClip(p0, p1, width, clip0, clip1, startT=0, endT=1):
    clipPlanes = GetRectoutClipPlanes(p0, p1, width, startT, endT)
    c0 = ClipAttributes()
    c0.quality = c0.Accurate  # Fast, Accurate
    c0.funcType = c0.Plane  # Plane, Sphere
    c0.plane1Status = 1
    c0.plane2Status = 1
    c0.plane3Status = 0
    c0.plane1Origin = clipPlanes[0][0]
    c0.plane2Origin = clipPlanes[1][0]
    c0.plane1Normal = clipPlanes[0][1]
    c0.plane2Normal = clipPlanes[1][1]
    c0.planeInverse = 1
    SetOperatorOptions(c0, clip0, 0)

    c1 = ClipAttributes()
    c1.quality = c1.Accurate  # Fast, Accurate
    c1.funcType = c1.Plane  # Plane, Sphere
    c1.plane1Status = 1
    c1.plane2Status = 1
    c1.plane3Status = 0 
    c1.plane1Origin = clipPlanes[2][0]
    c1.plane2Origin = clipPlanes[3][0]
    c1.plane1Normal = clipPlanes[2][1]
    c1.plane2Normal = clipPlanes[3][1]
    c1.planeInverse = 1
    SetOperatorOptions(c1, clip1, 0)

def Rectout(p0, p1, width, nbins):
    # Save the activeWindow
    activeWin = GetGlobalAttributes().activeWindow + 1

    # Create a copy of the current window
    CloneWindow()

    # Delete all but the first active plot
    pl = GetPlotList()
    activePlots = [1] * pl.GetNumPlots()
    plotvar = ""
    plotDB = ""
    for i in range(pl.GetNumPlots()):
        if pl.GetPlots(i).activeFlag == 1:
            plotvar = pl.GetPlots(i).plotVar
            plotDB = pl.GetPlots(i).databaseName
            activePlots[i] = 0
            break;
    if plotvar == "":
        DeleteWindow()
        SetActiveWindow(activeWin)
        raise "No active plots"
    if len(activePlots) > 1:
        SetActivePlots(activePlots)
        DeleteActivePlots()
    SetActivePlots(0)

    # Get the mesh var for the current plot
    md = GetMetaData(plotDB)
    for i in range(md.GetNumScalars()):
        if md.GetScalars(i).name == plotvar:
            varmesh = md.GetScalars(i).meshName 
    if varmesh == "":
        DeleteWindow()
        SetActiveWindow(activeWin)
        raise "Can't Rectout on non-scalar variable type"

    DefineScalarExpression("current_area", "area(%s)" % varmesh)
    DefineScalarExpression("original_area","area(%s)" % varmesh)
    DefineScalarExpression("scaled_var", "%s * current_area / original_area" % plotvar)

    # Add additional operators needed for the Rectout
    AddOperator("Clip")
    AddOperator("Clip")
    if AddOperator("DeferExpression") == 0:
        DeleteWindow()
        SetActiveWindow(activeWin)
        raise """
This script depends on the DeferExpression operator being loaded. You should 
make sure it is loaded in the Plugin Manager window and save your settings before
running this script again"""

    de = DeferExpressionAttributes(0)
    de.exprs = ("current_area", "scaled_var")
    SetOperatorOptions(de)

    # Get the indices of the last 2 Clip operators
    opNames = GetPlotList().GetPlots(0).operatorNames
    clipOps = []
    for i in range(len(opNames)):
        index = len(opNames)-1-i
        if opNames[index] == "Clip":
            clipOps.append(index)

    # Switch the active variable to the scaled_var
    ChangeActivePlotsVar("scaled_var")

    p0p1_len = vec_magnitude(vec_add(p1, vec_mult_float(p0, -1)))

    # Now, iterate over the bins and isolate each one in turn.
    binincr = 1. / float(nbins)
    curve = open("rectout.curve", "wt")
    curve.write("# rectout\n")
    for bin in range(nbins):
        start = float(bin) * binincr
        end = float(bin+1) * binincr
        RectClip(p0, p1, width, clipOps[-1], clipOps[-2], start, end)
        DrawPlots()

        # Distance along the segment p0 p1
        x = p0p1_len * (start + end) / 2.

        # Sum up the values in the area
        Query("Variable Sum")
        y = GetQueryOutputValue()

        # Save the values into the curve file
        curve.write("%g %g\n" % (x,y))
    curve.close()

    # Delete the temp window
    DeleteWindow()

    # Add a new curve plot window
    AddWindow()
    DeleteAllPlots()
    OpenDatabase("rectout.curve")
    AddPlot("Curve", "rectout")
    DrawPlots()

    # Go back to the original window
    SetActiveWindow(activeWin)

#
# Try an example.
#
def Example():
    OpenDatabase("~/data/rect2d.silo")
    AddPlot("Mesh", "curvmesh2d")
    AddPlot("Pseudocolor", "d")
    DrawPlots()
    Rectout((0.1666, 0.8333, 0), (0.8, 0.2, 0.), 0.25, 20)