Python Filters Development

Python filters Intro

VisIt does not currently have a way to directly manipulate datasets via python.

So, why do we care?

This type of functionality provides many potential benefits:

  • Runtime modification/prototyping of filters
  • Reduce development time for special purpose/one-off filters that are not of general use.
  • Path to powerful custom post processing

Python filters can play a useful role in following VisIt Realms:

  • Databases
  • Expressions
  • Operators
  • Queries

Development Tasks

Task Details Status
Embedded Python Interpreter & Python Filter Environment Create a simple embedded python interpreter in C++ with the ability to wrap/unwrap C++ VTK objects to/from python. Done
Create cmake option to enable python filters. We can't use them w/ static build, etc. Done
vtk-5.0.0d Prepare vtk distib with python wrappers. Done
build_visit changes Changed build_visit to build Python before VTK & enable VTK Python wrappers. Done
Python Expression Filter Add a basic python expression filter. Done
Python Query Filter Add a basic python query filter. Done
Python Source Editor Widget Create a simple python editor widget w/ syntax highlighting. (In progress)
Contract Module Create suitable python module for contract inspection & modification. (In progress)
Parallel Communication Module Create suitable python module for wrapping necessary mpi functions. See Mpicom Python Module for details. (In progress)
Python Operator Filter Add a basic python operator filter. Future Functionality
Use mpi4py Transition to use mpi4py for more functionality instead of small custom mpicom module. This will require some enhancements to mpi4py. Future Functionality

Development Notes

First Steps: VTK C++/Python Handshaking

VTK Provides Python wrappers for all of its classes, with these wrappers we can manipulate VTK objects across C++ / Python with a little effort.

C++ to Python

Python wrappers for VTK objects provide a constructor that optionally accepts a string for the address of a VTK object. This allows us to get the pointer from a C++ VTK object and pass it to an embedded python interpreter.

Example: 

// Create a python interpreter
PythonInterpreter pyi(argc,argv);
//  Setup / Manipulate your data set (or other VTK object)
vtkDataSet *ds = vtkDataSet::New();

// ...    

// Get a string representation of the pointer to the VTK object.
ostringstream oss;
oss << (void*) ds;
string addy_str = oss.str();
// Remove the "0x" prefix
addy_str = addy_str.substr(2);
// Feed the following into the python interpreter:
string py_script = "ds = vtk.vtkDataSet('" + addy_str + "')\n";
pyi.RunScript(py_script);
// "ds" is now available in the python interpreter

Python to C++

In python, you can obtain the address of any VTK object using the GetAddressAsString() method. It takes a string that specifies the vtk object name

Example:

In the Embedded Python Interpreter: 

# Setup / Manipulate your dataset (or other VTK Object)
ds = vtk.vtkUnstructuredGrid()

# ...

# Obtain the address string
ds_addy_str = ds.GetAddressAsString('vtkDataSet')
# Remove the prefix "Addy=" and convert it to a python integer
# The string is a hex address, so we need to pass the base to the 
# int conversion function
ds_addy = int(ds_addy_str[5:],16)

In C++ 

// Get access to the main python module
PyObject *main_mod = PyImport_AddModule((char*)"__main__");
// Get access to global dictionary
PyObject *global_d = PyModule_GetDict(main_mod);
// Get python object representing the address we obtained from the VTK Python object
PyObject *py_addy_int = PyDict_GetItemString(global_d, "ds_addy");
// Get the pointer value
long ds_addy = PyInt_AS_LONG(addy_int);
// Cast this pointer to our desired C++ instance 
vtkDataSet *ds = (vtkDataSet*)((void*)ds_addy);
// You know have access to the python created object in C++

Memory Magic

What happens to objects created in python when the interpreter is shutdown? Luckily we can increment the reference counter in C++ and when the interpreter is closed the object remains viable: 

vtkDataSet *ds = (vtkDataSet*)((void*)ds_addy);
// You now have access to the python created object in C++
ds->Register(NULL);
pyi.Shutdown();
// ds still valid

If you shutdown the interpreter before incrementing the reference count, do not expect a happy ending: 

vtkDataSet *ds = (vtkDataSet*)((void*)ds_addy);
// You now have access to the python created object in C++
pyi.Shutdown();
ds->Register(NULL);
// ds is invalid, expect a core dump

VisIt's version of VTK and VTK Python Wrappers

The VisIt team maintains its own version of VTK, a derivative of VTK 5.0.0. Doing so has allowed us to quickly incorporate critical bug fixes, exclude components not used by VisIt, and has shielded us from disruptive API changes. Please see VisIt VTK Modifications for a detailed list of what changes were made.

Prototype: Simple Single Input Python Expression

A prototype branch exists that provides the required changes to build_visit, a beta vtk-5.0.0d.tar.gz, and the implementation a simple python script expression. If you would like to take it for a spin please do so. Comments and suggestions would be very helpful.


Obtaining and setting up the branch

  • Checkout the branch
  co_branch 0804_python_filters src/ -user cyrush
  • Obtain vtk-5.0.0d.test.tar.gz
  svn export svn+ssh://myusername@portal-auth.nersc.gov/project/projectdirs/visit/svn/visit/branches/cyrush/0903_python_filters/third_party/vtk-5.0.0d.test.tar.gz
  • Run build_visit to build third party libraries - including the vtk python wrappers (make sure vtk-5.0.0d.test.tar.gz is in the working dir)
 src/svn_bin/build_visit --console --no-visit --no-thirdparty --cmake --mesa --python --vtk 
 (You should be able to use an existing Qt4 install)
  • Build the branch using the generated host.conf (sub in your qt build).

Python Script Expression Example

  • Launch VisIt
  • Open up rect2d.silo
  • Create an expression expr = py(quadmesh2d,"load_filter('simple.py')")
  • Create a file simple.py with the following source:
from math import sin, pi
def derive_var(ds_in):
    # ds_in is a vtk dataset, we want to create and return a new vtkDataArray
    ds_bounds = ds_in.GetBounds()
    # create a simple sine wave pattern
    x_ext = ds_bounds[1] - ds_bounds[0]
    y_ext = ds_bounds[3] - ds_bounds[2]
    
    ncells = ds_in.GetNumberOfCells()
    res = vtk.vtkFloatArray()
    res.SetNumberOfComponents(1)
    res.SetNumberOfTuples(ncells)
    for i in xrange(ncells):
        cell = ds_in.GetCell(i)
        bounds = cell.GetBounds()
        xv = bounds[0] + bounds[1] / 2.0
        yv = bounds[2] + bounds[3] / 2.0
        res.SetTuple1(i,.25 * (sin(xv*3*pi/x_ext) + sin(yv * 3*pi / y_ext)))
    return res

Note: You may also type the source directly into the expression window enclosed in quotes: 

  expr = py(quadmesh2d,"
  ... python ....
  ")

Create a Pseudocolor plot of "expr" and you should see the following:

Python filters simple example.png


Expression Prototype Example

# this expression works well for rect2d & rect3d.
from math import sin, pi
class MyExpression(SimplePythonExpression):
    def __init__(self):
        SimplePythonExpression.__init__(self)
        self.output_is_point_var  = False
        self.output_dimension = 1
    def derive_variable(self,ds_in,domain_id):
        # ds_in is a vtk dataset, we want
        # to create and return a new vtkDataArray
        # that contains a simple sine wave pattern 
        print ds_in, domain_id
        ds_bounds = ds_in.GetBounds()
        x_ext = ds_bounds[1] - ds_bounds[0]
        y_ext = ds_bounds[3] - ds_bounds[2]
        z_ext = ds_bounds[5] - ds_bounds[4]
        ncells = ds_in.GetNumberOfCells()
        res = vtk.vtkFloatArray()
        res.SetNumberOfComponents(1)
        res.SetNumberOfTuples(ncells)
        for i in xrange(ncells):
            cell = ds_in.GetCell(i)
            bounds = cell.GetBounds()
            xv = bounds[0] + bounds[1] / 2.0
            yv = bounds[2] + bounds[3] / 2.0
            zv = bounds[4] + bounds[5] / 2.0
            val = sin(xv*3*pi/x_ext) + sin(yv * 3*pi / y_ext)
            if z_ext != 0:
                val+= sin(zv * 3*pi / z_ext)
            res.SetTuple1(i,.2 * val)
        return res

py_filter = MyExpression


Expression Prototype Example (with MPI)

# this expression works well for rect2d & rect3d or multi_rect2d & multi_rect3d.
from math import sin, pi
class MyExpression(SimplePythonExpression):
    def __init__(self):
        SimplePythonExpression.__init__(self)
        self.output_is_point_var  = False
        self.output_dimension = 1
    def execute(self,dsets,domain_ids):
        # get bounds for entire dataset
        bounds = self.__get_bounds(dsets)
        self.x_ext = bounds[1] - bounds[0]
        self.y_ext = bounds[3] - bounds[2]
        self.z_ext = bounds[5] - bounds[4]
        # call base class method to do heavy lifting
        return SimplePythonExpression.execute(self,dsets,domain_ids)
    def derive_variable(self,ds_in,domain_id):
        ncells = ds_in.GetNumberOfCells()
        res = vtk.vtkFloatArray()
        res.SetNumberOfComponents(1)
        res.SetNumberOfTuples(ncells)
        print ncells
        for i in xrange(ncells):
            cell = ds_in.GetCell(i)
            bounds = cell.GetBounds()
            xv = (bounds[1] - bounds[0]) / 2.0  + bounds[0]
            yv = (bounds[3] - bounds[2]) / 2.0 + bounds[2]
            zv = (bounds[5] - bounds[4]) / 2.0 + bounds[4]
            #print xv,yv,zv
            wv = 2.0
            wv = wv * 2.0 * pi
            val = sin(xv*wv/self.x_ext) + sin(yv * wv/self.y_ext)
            if self.z_ext != 0.0:
                val+= sin(zv * wv/self.z_ext)
            res.SetTuple1(i,.2 * val)
        return res
    def __get_bounds(self,dsets):
        bounds = [ ds.GetBounds() for ds in dsets] 
        lbounds = [ 1e80] * 3
        ubounds = [-1e80] * 3
        if len(bounds) > 0:
            lbounds[0] = min(b[0] for b in bounds)
            lbounds[1] = min(b[2] for b in bounds)
            lbounds[2] = min(b[4] for b in bounds)
            ubounds[0] = max(b[1] for b in bounds)
            ubounds[1] = max(b[3] for b in bounds)
            ubounds[2] = max(b[5] for b in bounds)
        lbounds = mpicom.min(lbounds)
        ubounds = mpicom.max(ubounds)
        return [lbounds[0],ubounds[0],lbounds[1],ubounds[1],lbounds[2],ubounds[2]]
            
py_filter = MyExpression

Query Prototype Example

import mpicom
class CellAvgQuery(PythonQuery):
    def __init__(self):
        PythonQuery.__init__(self)
    def pre_execute(self):
        # init vars used to compute the average
        self.var_name     = self.attributes["VariableNames"][0]
        self.total_ncells = 0
        self.total_sum    = 0.0
    def execute(self,ds_in,domain_id):
        # sum over cell data array passed to query args
        ncells = ds_in.GetNumberOfCells()
        self.total_ncells += ncells
        cell_data = ds_in.GetCellData().GetArray(self.var_name)
        for i in xrange(ncells):
            self.total_sum += cell_data.GetTuple1(i)
    def post_execute(self):
        # communicate total sum & total # of cells
        # calculate average and set results
        self.total_ncells = mpicom.sum(self.total_cells)
        self.total_sum    = mpicom.sum(self.total_sum)
        result = self.total_sum / ((float)self.total_ncells)
        self.set_result_text("The average value of %s = %f" % (self.var_name,result))
        self.set_result_value(result)

py_filter = CellAvgQuery

Contract Wrappers

To give filters first class influence over the avt pipeline, they need to be able to modify avtContract instances. This requires some kind of python glue for the following avt classes:

  • avtContract
  • avtDataRequest
  • avtSILSpecification (shouldn't this be renamed to avtSILRequest?)
  • avtDataSelection

Note: avtDataSelection is basically a skeleton parent class, relevant subclasses need to be included as well.

Retrieved from "http://www.visitusers.org/index.php?title=Python_Filters_Development&oldid=5847"