VisIt-tutorial-Advanced-PythonQuery

VisIt's Python Query filter runtime allows you to use Python to process the data sets at the end of your visualization pipeline.

Resources

Data (from Tutorial Data):

• data/tgvortex.mir

Example Files (from Tutorial Examples):

• VisIt Python Client Driver: examples/advanced/python_query/visit_ccl_centroid_graph.py
• Python Query Source: examples/advanced/python_query/visit_ccl_centroid_graph.vpq

Processing connected components results with a Python Query

This tutorial will show how to create a Python query to further analyze connected components identified by VisIt.

The example Python Query demonstrates the following:

• Calculating an estimate of the centroid of each connected component.
• Constructing a simple graph of relationships between centroids.
• Saving the resulting graphs as VTK unstructured meshes for visualization.

To run the example:

• Launch VisIt and open the Taylor-Green vortex data set (data/tgvortex.mir).
• Launch the CLI or open the Commands window.
• Run the following script:

Source("examples/advanced/python_query/visit_ccl_centroid_graph.py")
setup_centroid_graphs()

This creates the following output files:

• The ccl_centroids_*.vtk files contain point meshes representing the centroids.
• The ccl_centorid_graph_*.vtk files contain line meshes with edges that connect each centroid to its closest neighboring centroid.

• Add a new window and clear any copied plots.
• Open output ccl_centroids_*.vtk
  • Add a Pseudocolor plot of component_label
  • Change Pseudocolor attributes: Point type: Sphere, Point size: 20
• Click Draw

This shows you the centroids of the identified connected components.

• Uncheck Apply operators to all plots
• Open output file ccl_centorid_graph_0000.vtk
  • Add a Pseudocolor plot of component_label
  • Accept the database correlation prompt
  • Add a Tube operator (Geometry/Tube)

• Click Draw

This shows a graph of the connected component centroids with edges connecting each centroid to its closest neighboring centroid.

• Lock the view and time sliders between both windows
• Accept the database correlation prompt
• Explore the results for the three time steps.

Example script files:

• Client Driver: examples/advanced/python_query/visit_ccl_centroid_graph.py

#
# file: visit_ccl_centroid_graph.py
#
# Driver which demonstrates using the python query "visit_ccl_centroid_graph.vpq"
#
# Usage:
#  Launch VisIt and open the "data/tgvortex.mir" file from tutorial data.
#  From the CLI run:
#  >>>Source("examples/advanced/python_query/visit_ccl_centroid_graph.py")
#  >>>setup_centroid_graphs()
#
import os
from os.path import join as pjoin

script_path = os.path.split(os.path.abspath(__visit_source_file__))[0]

def calc_centroid_graph(ts):
    """Execute our python query"""
    TimeSliderSetState(ts)
    Query("Max", use_actual_data=1)
    ncomps = int(GetQueryOutputValue()) + 1
    PythonQuery(file=pjoin(script_path,"visit_ccl_centroid_graph.vpq"),
	      vars=["cell_label","cell_volume"],
                args=[ts,ncomps])

def setup_ccl_plot():
    """Create connected components plot with cell volumes."""
    DeleteAllPlots()
    ReOpenDatabase(WindowInformation().activeSource)
    DefineScalarExpression("cell_label","conn_components(mesh)")
    DefineScalarExpression("cell_volume","volume(mesh)")
    AddPlot("Pseudocolor", "cell_label")
    AddOperator("Isovolume")
    iso_atts = IsovolumeAttributes()
    iso_atts.lbound = 0.8
    iso_atts.variable = "velocity_magnitude"
    SetOperatorOptions(iso_atts)
    AddOperator("DeferExpression")
    dexpr_atts = DeferExpressionAttributes()
    dexpr_atts.exprs = ("cell_label", "cell_volume")
    SetOperatorOptions(dexpr_atts)
    DrawPlots()
    
def setup_centroid_graphs():
    """ Main driver"""
    setup_ccl_plot()
    for ts in range(TimeSliderGetNStates()):
        calc_centroid_graph(ts)

• Python Query: examples/advanced/python_query/visit_ccl_centroid_graph.vpq

#
# file: visit_ccl_centroid_graph.vpq
#
# Creates a graph of CCL centroids with edges connecting 
# each centroid to its closest neighboring centroid.
#
# This example was created to be extra descriptive, and
# as a result it is not optimized. 
#
#

import math

def distance(c1,c2):
    """Euclidean Distance"""
    dx = c1[0] - c2[0]
    dy = c1[1] - c2[1]
    dz = c1[2] - c2[2]
    return math.sqrt(dx*dx+dy*dy+dz*dz)

class CCLCentroidGraph(SimplePythonQuery):
    """Custom Python Query"""
    def __init__(self):
        SimplePythonQuery.__init__(self)
        self.name = "CCLCentroidGraph"
        self.description = "Creating CCL Centroid Graph"
        self.obase = "ccl"
    def pre_execute(self):
        """Called prior to chunk execution."""
        # time step and number of components are passed in
        # from the python client
        self.ts = self.arguments[0]
        self.num_comps = self.arguments[1]
        self.cent_x = [0.0] * self.num_comps
        self.cent_y = [0.0] * self.num_comps
        self.cent_z = [0.0] * self.num_comps
        self.vol    = [0.0] * self.num_comps
    def execute_chunk(self,ds_in,domain_id):
        """
        Called to process each chunk of a domain decomposed mesh.
        Each MPI task processes a subset of the total chunks.
        """
        # get the cells in this chunk
        ncells = ds_in.GetNumberOfCells()
        # get the cell ccl labels, and the cell volumes
        cell_labels = ds_in.GetCellData().GetArray("cell_label")
        cell_vols   = ds_in.GetCellData().GetArray("cell_volume")
        # approximate centroid by integrating volume
        # scaled by bounding box center (per ccl label)
        for i in xrange(ncells):
            cell   = ds_in.GetCell(i)
            bounds = cell.GetBounds()
            cell_vol = cell_vols.GetTuple1(i)
            cell_lbl = int(cell_labels.GetTuple1(i))
            # bounding box center
            cell_bbx = (bounds[1] + bounds[0])/2.0
            cell_bby = (bounds[3] + bounds[2])/2.0
            cell_bbz = (bounds[5] + bounds[4])/2.0
            # volume scaling for integral estimate of the centroid 
            self.cent_x[cell_lbl] += cell_bbx * cell_vol
            self.cent_y[cell_lbl] += cell_bby * cell_vol
            self.cent_z[cell_lbl] += cell_bbz * cell_vol
            self.vol[cell_lbl]    += cell_vol
    def post_execute(self):
        """
        Called on all MPI tasks, after all chunks have been processed.
        """
        # get data from other procs
        self.cent_x = mpicom.sum(self.cent_x)
        self.cent_y = mpicom.sum(self.cent_y)
        self.cent_z = mpicom.sum(self.cent_z)
        self.vol    = mpicom.sum(self.vol)
        # scale by volume to obtain the final centroid estimate
        centroids = []
        for i in xrange(self.num_comps):
            if self.vol[i] != 0:
                self.cent_x[i] = self.cent_x[i] / self.vol[i]
                self.cent_y[i] = self.cent_y[i] / self.vol[i]
                self.cent_z[i] = self.cent_z[i] / self.vol[i]
                centroids.append([self.cent_x[i],self.cent_y[i],self.cent_z[i]])
        if mpicom.rank() == 0:
            # write our results
            self.write_centroid_points(centroids)
            self.write_centroid_graph(centroids)
            rtxt  = "Centroid Graph Results written to "
            rtxt += "%s and %s" % (self.points_ofile(),self.graph_ofile())
            self.set_result_text(rtxt)
            self.set_result_value(1)
    def points_ofile(self):
        return "%s_centroids_%04d.vtk" % (self.obase,self.ts)
    def graph_ofile(self):
        return "%s_centroid_graph_%04d.vtk" % (self.obase,self.ts)
    def write_centroid_points(self,centroids):
        """
        Create a vtk unstructured grid holding the centroid points,
        including the component label and the total component volume.
        """
        n_cents       = len(centroids)
        res           = vtk.vtkUnstructuredGrid()
        res_pts       = vtk.vtkPoints()
        # holds the component id of the centroid
        res_cent_lbls = vtk.vtkFloatArray()
        res_cent_lbls.SetNumberOfComponents(1)
        res_cent_lbls.SetNumberOfTuples(n_cents)
        res_cent_lbls.SetName("component_label")
        # holds the aggregate component volume
        res_cent_vols = vtk.vtkFloatArray()
        res_cent_vols.SetNumberOfComponents(1)
        res_cent_vols.SetNumberOfTuples(n_cents)
        res_cent_vols.SetName("component_volume")
        for i in xrange(n_cents):
            res_pts.InsertNextPoint(centroids[i])
            vertex = vtk.vtkVertex()
            vertex.GetPointIds().SetId(0, i)
            res.InsertNextCell(vertex.GetCellType(),vertex.GetPointIds())
            res_cent_lbls.SetTuple1(i,i)
            res_cent_vols.SetTuple1(i,self.vol[i])
        # complete our data set
        res.SetPoints(res_pts)
        res.GetCellData().AddArray(res_cent_lbls)
        res.GetCellData().AddArray(res_cent_vols)
        # write out the result
        writer = vtk.vtkUnstructuredGridWriter()
        writer.SetFileTypeToASCII()
        writer.SetInputData(res)
        writer.SetFileName(self.points_ofile())
        writer.Write()
    def write_centroid_graph(self,centroids):
        """ 
        Create a vtk unstructured grid that holds an embedded 
        graph of centroids. Edges connect each centroid to its 
        closest neighboring centroid.
        """
        n_cents     = len(centroids)
        res         = vtk.vtkUnstructuredGrid()
        res_pts     = vtk.vtkPoints()
        res_cent_lbls = vtk.vtkFloatArray()
        res_cent_lbls.SetNumberOfComponents(1)
        res_cent_lbls.SetNumberOfTuples(n_cents)
        res_cent_lbls.SetName("component_label")
        res_cent_vols = vtk.vtkFloatArray()
        res_cent_vols.SetNumberOfComponents(1)
        res_cent_vols.SetNumberOfTuples(n_cents)
        res_cent_vols.SetName("component_volume")
        for i in xrange(n_cents):
            res_pts.InsertNextPoint(centroids[i])
            res_cent_lbls.SetTuple1(i,i)
            res_cent_vols.SetTuple1(i,self.vol[i])
        # find closest centroid pairs
        cent_links = [-1] * n_cents
        for i in xrange(n_cents):
            if cent_links[i] == -1:
                dist = 1e32
                cid  = -1
                for j in xrange(n_cents):
                    if j != i:
                        curr_dist = distance(centroids[i],centroids[j])
                        if curr_dist < dist:
                            dist = curr_dist
                            cid  = j
                cent_links[i] = cid
                cent_links[j] = -2
        # create edges
        for i in xrange(n_cents):
            if cent_links[i] >= 0:
                line = vtk.vtkLine()
                line.GetPointIds().SetId(0, i)
                line.GetPointIds().SetId(1, cent_links[i])
                res.InsertNextCell(line.GetCellType(),line.GetPointIds())
        # complete our data set
        res.SetPoints(res_pts)
        res.GetPointData().AddArray(res_cent_lbls)
        res.GetPointData().AddArray(res_cent_vols)
        # write out the result
        writer = vtk.vtkUnstructuredGridWriter()
        writer.SetFileTypeToASCII()
        writer.SetInputData(res)
        writer.SetFileName(self.graph_ofile())
        writer.Write()
 
py_filter = CCLCentroidGraph