VTK tools

Pygslib use VTK:

  • as data format and data converting tool
  • to plot in 3D
  • as a library with some basic computational geometry functions, for example to know if a point is inside a surface

Some of the functions in VTK were obtained or modified from Adamos Kyriakou at https://pyscience.wordpress.com/


In [1]:
import pygslib 
import numpy as np


Functions in vtktools


In [2]:
help(pygslib.vtktools)


Help on module pygslib.vtktools in pygslib:

NAME
    pygslib.vtktools - PyGSLIB vtktools, Module with tools to work with VTK.

FILE
    c:\og_python\pygslib\pygslib\vtktools.pyd

DESCRIPTION
    Copyright (C) 2015 Adrian Martinez Vargas 
    
    This program is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 3 of the License, or
    any later version.
       
    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.
       
    You should have received a copy of the GNU General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>

FUNCTIONS
    GetPointsInPolydata(...)
        p= GetPointsInPolydata(polydata)
        
        Returns the points coordinates in a polydata (e.j. wireframe) 
        as a numpy array. 
        
        
        Parameters
        ----------
        polydata : VTK polydata
        
            
        Returns
        -------
        p : 2D array of float
            Points coordinates in polydata object  
            
        See also
        --------
        SetPointsInPolydata, SavePolydata
    
    SavePolydata(...)
        SavePolydata(polydata, path)
        
        Save polydata in a VTK XML polydata file (ext vtp) 
        
        
        Parameters
        ----------
        polydata : VTK polydata
        path : str 
                Extension (*.vtp) will be added if not provided
    
    SetPointsInPolydata(...)
        SetPointsInPolydata(polydata)
        
        Set the points coordinates in a polydata (e.j. wireframe). 
        
        
        Parameters
        ----------
        polydata : VTK polydata
        points   : 2D array 
                New points coordinates to be set in the polydata
                points shape may be equal to [polydata.GetNumberOfPoints(),3] 
        
            
        See also
        --------
        GetPointsInPolydata, SavePolydata
    
    addLine(...)
        addLine(renderer, 
                 p1,
                 p2,  
                 color=(0.0, 0.0, 1.0))
        
        Adds a line into an existing VTK renderer 
        
        Parameters
        ----------
        renderer : VTK renderer
               renderer with vtk objects and properties 
        p1,p2: tuple with 3 float
               point location
        radius: radius of the point  
               Default 1.
        color: tuple with 3 float
                Default (0.0, 0.0, 0.0)
        
        Returns
        -------
        VtkRenderer : an Vtk renderer containing VTK polydata 
           
        See Also
        --------
        vtk_show, loadSTL
        
        Notes
        -----
        This Code was modified fron the original code by Adamos Kyriakou
        published in https://pyscience.wordpress.com/
    
    addPoint(...)
        addPoint(renderer, 
                 p, 
                 radius=1.0, 
                 color=(0.0, 0.0, 0.0))
        
        Adds a point into an existing VTK renderer 
        
        Parameters
        ----------
        renderer : VTK renderer
               renderer with vtk objects and properties 
        p: tuple with 3 float
               point location
        radius: radius of the point  
               Default 1.
        color: tuple with 3 float
                Default (0.0, 0.0, 0.0)
        
        Returns
        -------
        VtkRenderer : an Vtk renderer containing VTK polydata 
           
        See Also
        --------
        vtk_show, loadSTL
        
        Notes
        -----
        This Code was modified fron the original code by Adamos Kyriakou
        published in https://pyscience.wordpress.com/
    
    dmtable2wireframe(...)
        mywireframe = dmtable2wireframe(x, y, z, pid1,pid2,pid3) 
        
        Takes a wireframe defined by two table: 
        
        - x, y, z : This is the points table
        - pid1,pid2,pid3: This is another table defining triangles with 
          three point IDs. 
            
        Parameters
        ---------- 
        x,y,z          : 1D array of floats
            coordinates of the points to be tested        
        pid1,pid2,pid3 : 1D array of integers
            triangles defined by three existing point IDs
        str filename   : Str (Default None)
            file name and path. If provided the file wireframe is saved to 
            this file.   
        
        Returns
        -------
        surface : VTK polydata
               wireframe imported
        
        
        Notes
        -----
        ``pid`` is the row number in the point table.
        
        ``pid`` indices start at zero
        
        For example: 
        
        a point table
        
        | x | y | z |
        -------------
        | 0 | 0 | 0 |
        | 1 | 0 | 0 | 
        | 0 | 1 | 0 |
        
        a triangle table
        
        |pid1|pid2|pid3|
        -------------
        | 0  | 1  | 2  |
    
    getbounds(...)
        boundingbox(polydata)
        
        Returns the bounding box limits x,y,z minimum and maximum
        
        
        Parameters
        ----------
        polydata : VTK polydata
        
            
        Returns
        -------
        xmin,xmax, ymin,ymax, zmin,zmax : float
            The geometry bounding box
    
    grid2vtkfile(...)
        grid2vtkfile(path, x,y,z, data)
        
        save grid in vtk file if x,y,z are 1D it saves VtkRectilinearGrid,
        if are 3D it saves VtkStructuredGrid
        
        
        Parameters 
        ----------
        
        
            
        Returns
        -------
    
    loadSTL(...)
        loadSTL(filenameSTL)
        
        Load a STL wireframe file
        
        Parameters
        ----------
        filenameSTL : file path
        
        Returns
        -------
        polydata : VTK Polydata object 
               
           
        Notes
        -----
        Code by Adamos Kyriakou
        published in https://pyscience.wordpress.com/
    
    partialgrid2vtkfile(...)
        partialgrid2vtkfile(path, x,y,z, dx,dy,dz,var, varname)
        
        save unstructurw vtk grid of parent blocks into a file
        
        
        Parameters
        ----------
        
        
            
        Returns
        -------
    
    pointinsolid(...)
        inside = pygslib.vtktools.pointinsolid(surface, x, y, z)
        
        Find points inside a closed surface using vtkSelectEnclosedPoints
            
        
        Parameters
        ----------
        surface : VTK polydata
               this may work with any 3D object..., no necessarily triangles   
        x,y,z   : 1D array of floats
               coordinates of the points to be tested
          
        
        Returns
        -------
        inside : 1D array of integers 
            Indicator of point inclusion with values [0,1]  
            0 means that the point is not inside
            1 means that the point is inside
        
        See Also
        --------
        vtk_raycasting, pointquering
        
        Notes
        -----
        It is assumed that the surface is closed and manifold.  
        Note that if this check is not performed and the surface is not 
        closed, the results are undefined.
    
    pointquering(...)
        query = pygslib.vtktools.pointquering(surface, azm, dip, x, y, z, test)
        
        Find points inside, over and below surface/solid
            
        
        Parameters
        ----------
        surface : VTK polydata
               this may work with any 3D object..., no necessarily triangles   
        azm, dip: float
               rotation defining the direction we will use to test the points
               azm 0 will point north and dip positive meas downward direction
               (like surface drillholes)
        x,y,z   : 1D array of floats
               coordinates of the points to be tested
        test    : integer
               1 test inside closed surface. Here we use 
                 vtkOBBTree::InsideOrOutside. Closed surface are required
               2 test 'above' surface 
               3 test 'below' surface 
               4 test 'inside' surface (the surface can be open)
        
        Returns
        -------
        inside : 1D array of integers 
            Indicator of point inclusion with values [0,1]  
            0 means that the point is not inside, above or below surface
            1 means that the point is inside, above or below surface
        p1 : 1D array size 3
            rays to show where are pointing in the output
        
        
        See Also
        --------
        vtk_raycasting
        
        Notes
        -----
        The test 1 requires the surface to be close, to find points between
        two surfaces you can use test=4
        The test, 2,3 and 4 use raycasting and rays orientation are defined
        by azm and dip values. The results will depend on this direction,
        for example, to find the points between two open surfaces you may 
        use ray direction perpendicular to surfaces  
        If a surface is closed the points over and below a surface will 
        include points inside surface
    
    points2vtkfile(...)
        bpoints2vtkfile(path, x,y,z, data)
        
        Save points in vtk file
        
        
        Parameters
        ----------
        
        
            
        Returns
        -------
    
    polydata2renderer(...)
        polydata2renderer(polydata, 
                 color=None, 
                 opacity=None, 
                 background=None)
        
        Creates vtk renderer from vtk polydata
        
        Parameters
        ----------
        polydata : VTK polydata
        
        color: tuple with 3 floats (RGB)
            default color=(1.,0.,0.)
        opacity: float 
            default 1 (opaque)
        background: tuple with 3 floats (RGB)
            default background=(1.,1.,1.)
            
            
        Returns
        -------
        VtkRenderer : an Vtk renderer containing VTK polydata 
           
        See Also
        --------
        vtk_show, loadSTL
        
        Notes
        -----
        This Code was modified fron the original code by Adamos Kyriakou
        published in https://pyscience.wordpress.com/
    
    vtk_raycasting(...)
        vtk_raycasting(surface, pSource, pTarget)
        
        Intersects a line defined by two points with a polydata vtk object,
        for example a closed surface. 
        
        This function is to find intersection and inclusion of points 
        within, below or over a surface. This is handy to select points 
        or define blocks inside solids.    
        
        Parameters
        ----------
        surface : VTK polydata
               this may work with any 3D object..., no necessarily triangles   
        pSource, pTarget: tuples with 3 float
               point location defining a ray
        
        Returns
        -------
        intersect : integer, with values -1, 0, 1
            0 means that no intersection points were found
            1 means that some intersection points were found
           -1 means the ray source lies within the surface
        pointsIntersection : tuple of tuples with 3 float
            this is a tuple with coordinates tuples defining 
            intersection points
        
        pointsVTKIntersectionData : an Vtk polydata object 
            similar to pointsIntersection but in Vtk polydata format
        
        
        See Also
        --------
        vtk_show, loadSTL
        
        Notes
        -----
        This Code was modified fron the original code by Adamos Kyriakou
        published in https://pyscience.wordpress.com/
    
    vtk_show(...)
        vtk_show(renderer, 
                 width=400, 
                 height=300,
                 camera_position=None,
                 camera_focalpoint=None)
        
        Display a vtk renderer in Ipython Image
        
        Parameters
        ----------
        renderer : VTK renderer
               renderer with vtk objects and properties 
        width, height: float
               Size of the image, default (400,300)
        camera_position, camera_focalpoint: tuples with 3 float 
               camera position, camera focal point (ex. center of mass)
               default None. If not None the the camera will be overwrite
        Returns
        -------
        Image : an IPython display Image object
          
        See Also
        --------
        polydata2renderer, loadSTL
           
        Notes
        -----
        This Code was modified fron the original code by Adamos Kyriakou
        published in https://pyscience.wordpress.com/

DATA
    __test__ = {}


Load a cube defined in an stl file and plot it

STL is a popular mesh format included an many non-commercial and commercial software, example: Paraview, Datamine Studio, etc.


In [3]:
#load the cube 
mycube=pygslib.vtktools.loadSTL('../datasets/stl/cube.stl')

# see the information about this data... Note that it is an vtkPolyData
print mycube


vtkPolyData (00000000114BF7D0)
  Debug: Off
  Modified Time: 185
  Reference Count: 1
  Registered Events: (none)
  Information: 0000000007BD7FC0
  Data Released: False
  Global Release Data: Off
  UpdateTime: 186
  Field Data:
    Debug: Off
    Modified Time: 146
    Reference Count: 1
    Registered Events: (none)
    Number Of Arrays: 0
    Number Of Components: 0
    Number Of Tuples: 0
  Number Of Points: 8
  Number Of Cells: 12
  Cell Data:
    Debug: Off
    Modified Time: 154
    Reference Count: 1
    Registered Events: 
      Registered Observers:
        vtkObserver (0000000007FC9AE0)
          Event: 33
          EventName: ModifiedEvent
          Command: 0000000002BA42A0
          Priority: 0
          Tag: 1
    Number Of Arrays: 0
    Number Of Components: 0
    Number Of Tuples: 0
    Copy Tuple Flags: ( 1 1 1 1 1 0 1 1 )
    Interpolate Flags: ( 1 1 1 1 1 0 0 1 )
    Pass Through Flags: ( 1 1 1 1 1 1 1 1 )
    Scalars: (none)
    Vectors: (none)
    Normals: (none)
    TCoords: (none)
    Tensors: (none)
    GlobalIds: (none)
    PedigreeIds: (none)
    EdgeFlag: (none)
  Point Data:
    Debug: Off
    Modified Time: 156
    Reference Count: 1
    Registered Events: 
      Registered Observers:
        vtkObserver (0000000007FC98D0)
          Event: 33
          EventName: ModifiedEvent
          Command: 0000000002BA42A0
          Priority: 0
          Tag: 1
    Number Of Arrays: 0
    Number Of Components: 0
    Number Of Tuples: 0
    Copy Tuple Flags: ( 1 1 1 1 1 0 1 1 )
    Interpolate Flags: ( 1 1 1 1 1 0 0 1 )
    Pass Through Flags: ( 1 1 1 1 1 1 1 1 )
    Scalars: (none)
    Vectors: (none)
    Normals: (none)
    TCoords: (none)
    Tensors: (none)
    GlobalIds: (none)
    PedigreeIds: (none)
    EdgeFlag: (none)
  Bounds: 
    Xmin,Xmax: (-5, 5)
    Ymin,Ymax: (-5, 5)
    Zmin,Zmax: (-5, 5)
  Compute Time: 200
  Number Of Points: 8
  Point Coordinates: 0000000010779F10
  Locator: 0000000000000000
  Number Of Vertices: 0
  Number Of Lines: 0
  Number Of Polygons: 12
  Number Of Triangle Strips: 0
  Number Of Pieces: 1
  Piece: 0
  Ghost Level: 0



In [4]:
# Create a VTK render containing a surface (mycube)
renderer = pygslib.vtktools.polydata2renderer(mycube, color=(1,0,0), opacity=0.50, background=(1,1,1))
# Now we plot the render
pygslib.vtktools.vtk_show(renderer, camera_position=(-20,20,20), camera_focalpoint=(0,0,0))


Out[4]:

Ray casting to find intersections of a lines with the cube

This is basically how we plan to find points inside solid and to define blocks inside solid


In [5]:
# we have a line, for example a block model row 
# defined by two points or an infinite line passing trough a dillhole sample 
pSource = [-50.0, 0.0, 0.0]
pTarget = [50.0, 0.0, 0.0]

# now we want to see how this looks like
pygslib.vtktools.addLine(renderer,pSource, pTarget, color=(0, 1, 0))
pygslib.vtktools.vtk_show(renderer) # the camera position was already defined


Out[5]:

In [6]:
# now we find the point coordinates of the intersections 
intersect, points, pointsVTK= pygslib.vtktools.vtk_raycasting(mycube, pSource, pTarget)

print "the line intersects? ", intersect==1
print "the line is over the surface?", intersect==-1

# list of coordinates of the points intersecting 
print points


the line intersects?  True
the line is over the surface? False
[(-5.0, 0.0, 0.0), (5.0, 0.0, 0.0)]

In [7]:
#Now we plot the intersecting points

# To do this we add the points to the renderer 
for p in points: 
    pygslib.vtktools.addPoint(renderer, p, radius=0.5, color=(0.0, 0.0, 1.0))

pygslib.vtktools.vtk_show(renderer)


Out[7]:

Test line on surface


In [8]:
# we have a line, for example a block model row 
# defined by two points or an infinite line passing trough a dillhole sample 
pSource = [-50.0, 5.01, 0]
pTarget = [50.0, 5.01, 0]

# now we find the point coordinates of the intersections 
intersect, points, pointsVTK= pygslib.vtktools.vtk_raycasting(mycube, pSource, pTarget)

print "the line intersects? ", intersect==1
print "the line is over the surface?", intersect==-1

# list of coordinates of the points intersecting 
print points

# now we want to see how this looks like
pygslib.vtktools.addLine(renderer,pSource, pTarget, color=(0, 1, 0))

for p in points: 
    pygslib.vtktools.addPoint(renderer, p, radius=0.5, color=(0.0, 0.0, 1.0))

pygslib.vtktools.vtk_show(renderer) # the camera position was already defined 

# note that there is a tolerance of about 0.01


the line intersects?  True
the line is over the surface? False
[(-5.0, 5.010000228881836, 0.0), (5.0, 5.010000228881836, 0.0)]
Out[8]:

Finding points


In [9]:
#using same cube but generation arbitrary random points
x = np.random.uniform(-10,10,150)
y = np.random.uniform(-10,10,150)
z = np.random.uniform(-10,10,150)

Find points inside a solid


In [10]:
# selecting all inside the solid
# This two methods are equivelent but test=4 also works with open surfaces  
inside,p=pygslib.vtktools.pointquering(mycube, azm=0, dip=0, x=x, y=y, z=z, test=1)
inside1,p=pygslib.vtktools.pointquering(mycube, azm=0, dip=0, x=x, y=y, z=z, test=4)
err=inside==inside1
#print inside, tuple(p)
print x[~err]
print y[~err]
print z[~err]


[]
[]
[]

In [11]:
# here we prepare to plot the solid, the x,y,z indicator and we also 
# plot the line (direction) used to ray trace

# convert the data in the STL file into a renderer and then we plot it
renderer = pygslib.vtktools.polydata2renderer(mycube, color=(1,0,0), opacity=0.70, background=(1,1,1))
# add indicator (r->x, g->y, b->z)
pygslib.vtktools.addLine(renderer,[-10,-10,-10], [-7,-10,-10], color=(1, 0, 0))
pygslib.vtktools.addLine(renderer,[-10,-10,-10], [-10,-7,-10], color=(0, 1, 0))
pygslib.vtktools.addLine(renderer,[-10,-10,-10], [-10,-10,-7], color=(0, 0, 1))

# add ray to see where we are pointing
pygslib.vtktools.addLine(renderer, (0.,0.,0.), tuple(p), color=(0, 0, 0))


Out[11]:
(vtkOpenGLRenderer)00000000114F7780

In [12]:
# here we plot the points selected and non-selected in different color and size
# add the points selected
for i in range(len(inside)):
    p=[x[i],y[i],z[i]]
    
    if inside[i]!=0:
        #inside
        pygslib.vtktools.addPoint(renderer, p, radius=0.5, color=(0.0, 0.0, 1.0))
    else:
        pygslib.vtktools.addPoint(renderer, p, radius=0.2, color=(0.0, 1.0, 0.0))

        
#lets rotate a bit this
pygslib.vtktools.vtk_show(renderer, camera_position=(0,0,50), camera_focalpoint=(0,0,0))


Out[12]:

Find points over a surface


In [13]:
# selecting all over a solid (test = 2) 
inside,p=pygslib.vtktools.pointquering(mycube, azm=0, dip=0, x=x, y=y, z=z, test=2)

# here we prepare to plot the solid, the x,y,z indicator and we also 
# plot the line (direction) used to ray trace

# convert the data in the STL file into a renderer and then we plot it
renderer = pygslib.vtktools.polydata2renderer(mycube, color=(1,0,0), opacity=0.70, background=(1,1,1))
# add indicator (r->x, g->y, b->z)
pygslib.vtktools.addLine(renderer,[-10,-10,-10], [-7,-10,-10], color=(1, 0, 0))
pygslib.vtktools.addLine(renderer,[-10,-10,-10], [-10,-7,-10], color=(0, 1, 0))
pygslib.vtktools.addLine(renderer,[-10,-10,-10], [-10,-10,-7], color=(0, 0, 1))

# add ray to see where we are pointing
pygslib.vtktools.addLine(renderer, (0.,0.,0.), tuple(-p), color=(0, 0, 0))


Out[13]:
(vtkOpenGLRenderer)000000001150E728

In [14]:
# here we plot the points selected and non-selected in different color and size
# add the points selected
for i in range(len(inside)):
    p=[x[i],y[i],z[i]]
    
    if inside[i]!=0:
        #inside
        pygslib.vtktools.addPoint(renderer, p, radius=0.5, color=(0.0, 0.0, 1.0))
    else:
        pygslib.vtktools.addPoint(renderer, p, radius=0.2, color=(0.0, 1.0, 0.0))

        
#lets rotate a bit this
pygslib.vtktools.vtk_show(renderer, camera_position=(0,0,50), camera_focalpoint=(0,0,0))


Out[14]:

Find points below a surface


In [15]:
# selecting all over a solid (test = 2) 
inside,p=pygslib.vtktools.pointquering(mycube, azm=0, dip=0, x=x, y=y, z=z, test=3)

# here we prepare to plot the solid, the x,y,z indicator and we also 
# plot the line (direction) used to ray trace

# convert the data in the STL file into a renderer and then we plot it
renderer = pygslib.vtktools.polydata2renderer(mycube, color=(1,0,0), opacity=0.70, background=(1,1,1))
# add indicator (r->x, g->y, b->z)
pygslib.vtktools.addLine(renderer,[-10,-10,-10], [-7,-10,-10], color=(1, 0, 0))
pygslib.vtktools.addLine(renderer,[-10,-10,-10], [-10,-7,-10], color=(0, 1, 0))
pygslib.vtktools.addLine(renderer,[-10,-10,-10], [-10,-10,-7], color=(0, 0, 1))

# add ray to see where we are pointing
pygslib.vtktools.addLine(renderer, (0.,0.,0.), tuple(p), color=(0, 0, 0))


Out[15]:
(vtkOpenGLRenderer)0000000011DDB6D0

In [16]:
# here we plot the points selected and non-selected in different color and size
# add the points selected
for i in range(len(inside)):
    p=[x[i],y[i],z[i]]
    
    if inside[i]!=0:
        #inside
        pygslib.vtktools.addPoint(renderer, p, radius=0.5, color=(0.0, 0.0, 1.0))
    else:
        pygslib.vtktools.addPoint(renderer, p, radius=0.2, color=(0.0, 1.0, 0.0))

        
#lets rotate a bit this
pygslib.vtktools.vtk_show(renderer, camera_position=(0,0,50), camera_focalpoint=(0,0,0))


Out[16]:

Export points to a VTK file


In [17]:
data = {'inside': inside}

pygslib.vtktools.points2vtkfile('points', x,y,z, data)

The results can be ploted in an external viewer, for example mayavi or paraview:


In [ ]: