In this post I will demonstrate how to use VTK to perform ray-casting, i.e., intersecting lines/rays with surface meshes, and retrieving the coordinates of those intersection points. This will include loading such meshes from STL files, performing the intersections, as well as visualizing the mesh, lines, and points in VTK. In addition, I will present my very own
pycaster package which cleanly wraps the VTK parts and allows for cleaner code and added functionality.
Today I’ll be talking about ray-casting which, to quote Wikipedia, is “the use of ray-surface intersection tests to solve a variety of problems in computer graphics and computational geometry”.
A pivotal difference between ray-casting and ‘ray-tracing’ is that the former only ‘casts’ a single ray, tests for its intersection with objects, and stops there. Ray-tracing on the other hand, is more physically accurate as it treats the rays with physics laws on reflection, refraction, attenuation, etc to ‘trace’, i.e., follow, that ray and its derivative rays.
However, I should make clear that ray-casting is the natural precursor to ray-tracing as it tells us what part of which object the ray intersects with and provides all necessary information to cast subsequent rays.
I will first demonstrate the approach on a very simple model of a hollow sphere which was designed in Rhino3D and exported to STL. Upon establishing the approach, I will talk about my very own
pycaster package which I wrote to wrap the VTK part of the code, simplify the whole process, and provide some additional functionality. Subsequently, I will repeat the ray-casting process with
pycaster on the model of the human skull we extracted and saved as an STL file in the previous post.
As I keep writing new posts I will try my darndest to build upon previously presented snippets in the interest of code-reuse, consistency, and brevity. I will be utilizing those snippets in the form of ‘helper-functions’, briefly documenting their function and referencing past posts where they were presented in more detail.
The following helper-functions will be used:
vtk_show(renderer, width=400, height=300): This function allows me to pass a
vtkRendererobject and get an IPython Notebook compatible image output of that render with given dimensions. This code was presented in this past post about VTK integration with an IPython Notebook.
addPoint(renderer, p, radius=1.0, color=[0.0, 0.0, 0.0]): This function allows me to pass the aforementioned
vtkRendererobject, which ‘holds’ the entire scene, and to add a ‘point’ to that scene (prior to its being rendered using
vtk_show). Upon inspecting the code you will see that we’re in fact rendering a sphere with its center coordinates stored within a
pparameter. The sphere’s
1.0, while we can also control its
colorwhich defaults to the RGB value for black. This code was presented in this past post about VTK integration with an IPython Notebook, where I showed how to render a red sphere.
addLine(renderer, p1, p2, color=[0.0, 0.0, 1.0]): This function is very similar to
addPointbut instead adds a line with
p2coordinates using the
vtkLineSourceclass. Its code is truly similar to
addPointso just take a look at this post’s IPython Notebook and this past post to understand the mechanics.
loadSTL(filenameSTL): This function uses the
vtkSTLReaderclass to load the contents of an
filenameSTLand return a
vtkPolyDataobject containing the mesh. While I’ve demonstrated how to write an STL file in the previous post, I haven’t exactly presented this function so here’s its code:
def loadSTL(filenameSTL): readerSTL = vtk.vtkSTLReader() readerSTL.SetFileName(filenameSTL) # 'update' the reader i.e. read the .stl file readerSTL.Update() polydata = readerSTL.GetOutput() # If there are no points in 'vtkPolyData' something went wrong if polydata.GetNumberOfPoints() == 0: raise ValueError( "No point data could be loaded from '" + filenameSTL) return None return polydata
Similarly to how one writes an STL file (shown in this previous post), we simply create a new
vtkSTLReader object, pass the filename of the
.stl file through the
SetFileName method, and call
Update to read in that file.
The mesh is then stored internally in a
vtkPolyData object, a pointer to which we retrieve through
GetOutput and store under
polydata which we eventually return. A point of interest in the above code lies in the usage of the
GetNumberOfPoints method of the
vtkPolyData class through which we check if the
vtkSTLReader object managed to read in anything. As I’ve said time and time again, VTK is not big on throwing exceptions or giving warnings. Thus, one should always check that the operation was indeed successful.
Ray-Casting with Python & VTK
In this example I will show the pure Python+VTK code so you can understand the entire process, while the simpler
pycaster version will follow . You can find the corresponding notebook here, while the STL model of the hollow-sphere can be download here.
Loading and rendering the mesh
Initially we load and render the mesh-surface of the hollow sphere using the
vtk_show helper-functions presented prior:
mesh = loadSTL("sphereHollow.stl") mapper = vtk.vtkPolyDataMapper() mapper.SetInput(mesh) actor = vtk.vtkActor() actor.SetMapper(mapper) actor.GetProperty().SetOpacity(0.25) renderer = vtk.vtkRenderer() renderer.AddActor(actor) renderer.SetBackground(1.0, 1.0, 1.0) vtk_show(renderer)
If you’re still having trouble getting this code then do read the past post about VTK integration with an IPython Notebook. The only novelty here is the
GetProperty method of the
vtkActor class. This method returns a pointer to this actor’s
vtkProperty object which we can then use to set a whole bunch of visual properties such as the color, visibility of edges, lighting, shading, etc. In this case, we use it to set this actor’s ‘opacity’ through the
SetOpacity method to
25% so we can see ‘inside’ the surface. Lastly, I want to emphasize that the
vtkRenderer is stored under the
renderer variable to which we’ll keep adding actors as we proceed while
mesh holds a pointer to the
vtkPolyData loaded from the STL file. The above code yields the below figure.
Preparing and rendering the ray
This hollow sphere model is centered around
(0,0,0), exhibits an outer radius of
25 mm, and an inner radius of
20 mm, resulting in a shell thickness of
5 mm. Therefore, for the purposes of this example, we set the source and target points of the ray to
50 mm on each side of the the
x axis as such:
pSource = [-50.0, 0.0, 0.0] pTarget = [50.0, 0.0, 0.0]
Note that there’s nothing special about the above ‘points’. They’re simply two
listobjects with three
To make the whole thing more visual we use the
addLine helper-functions to add two actors representing the ray’s points and an actor representing the line/ray to the
addPoint(renderer, pSource, color=[1.0, 0.0, 0.0]) addPoint(renderer, pTarget, color=[0.0, 1.0, 0.0]) addLine(renderer, pSource, pTarget) vtk_show(renderer)
As you can see the ‘source’ point will be rendered as red, while the ‘target’ is green. After adding the three actors, we simply use
vtk_show again and get the following figure.
Testing for intersection
Now here comes the fun part! At this point we have the surface-mesh under the
mesh variable and we’ve configured the ray, defined by the coordinates in
pTarget, we want to test for intersection with. The heart of this approach, however, lies in the usage of
vtkOBBTree which generates a oriented bounding-box ‘tree’ for the given geometry (the mesh in our case).
Firstly we need to prepare the OBB tree as follows:
obbTree = vtk.vtkOBBTree() obbTree.SetDataSet(mesh) obbTree.BuildLocator()
That’s it! We now have a world-class intersection tester at our disposal. We create a new
obbTree, set our
mesh as its dataset through
SetDataSet, and call
BuildLocator which will create the OBB tree, and allow for super-fast intersection testing :).
Before we proceed I want you to check the docs for the
vtkOBBTree. The method we’ll be using is
IntersectWithLine which as you can see sports 4 overloaded versions. The signature of the one we’ll be using is the following:
int vtkOBBTree::IntersectWithLine(const double a0, const double a1, vtkPoints * points, vtkIdList * cellIds)
which provides the following docstring: “Take the passed line segment and intersect it with the data set. This method assumes that the data set is a vtkPolyData that describes a closed surface, and the intersection points that are returned in ‘points’ alternate between entrance points and exit points. The return value of the function is 0 if no intersections were found, –1 if point ‘a0’ lies inside the closed surface, or +1 if point ‘a0’ lies outside the closed surface. Either ‘points’ or ‘cellIds’ can be set to NULL if you don’t want to receive that information.”
As in the case of ray-casting we don’t care about the cells with which the ray will intersect but do want information on the intersection points, we call the aforementioned method as such:
pointsVTKintersection = vtk.vtkPoints() code = obbTree.IntersectWithLine(pSource, pTarget, pointsVTKintersection, None)
As you can see we create a new
vtkPoints object, as required by the
IntersectWithLine signature, which we call
pointsVTKintersection. We then pass the
pTarget coordinates of the ray along with
pointsVTKintersection where the intersection points will be stored. Lastly, we’re using
cellIds which will be interpreted as
NULL as we don’t want that information.
codevariable which stores the return value of
IntersectWithLine? As you saw in this method’s docstring, a particular ‘code’ is returned depending on the result of the intersection test. It’d be a very good idea to always ‘check’ that value and ensure that intersection points were indeed found.
Extracting the intersection point coordinates
At this point we have acquired the intersection points which are stored under
pointsVTKintersection, an object of
vtkPoints type. However, we still need to ‘convert’ those
vtkPoints to something that can be manipulated in Python. This next part of the code is important:
pointsVTKIntersectionData = pointsVTKintersection.GetData() noPointsVTKIntersection = pointsVTKIntersectionData.GetNumberOfTuples() pointsIntersection =  for idx in range(noPointsVTKIntersection): _tup = pointsVTKIntersectionData.GetTuple3(idx) pointsIntersection.append(_tup)
As you can see, we first need to ‘extract’ the actual
pointsVTKintersection using the
GetData method, which we then assign to
pointsVTKIntersectionData (lengthy name, I know). We then get the number of intersection points found through
pointsVTKIntersectionData.GetNumberOfTuples() as every point will be stored as a 3-value tuple of the 3D coordinates. We create an empty
pointsIntersection which will store all those tuples. Subsequently, we loop through those tuples stored in
pointsVTKIntersectionData, acquiring each of them through the
GetTuple3 method, appending each tuple in the
pointsIntersection list. A simple
[(-24.91461181640625, 0.0, 0.0), (-19.914684295654297, 0.0, 0.0), (19.914684295654297, 0.0, 0.0), (24.91461181640625, 0.0, 0.0)]
Finally, using the
vtk_show helper functions we loop through these points and render them (in blue) as such:
for p in pointsIntersection: addPoint(renderer, p, color=[0.0, 0.0, 1.0]) vtk_show(renderer)
yielding the figure below.
pycaster in order to wrap and simplify the VTK functionality you saw above. The
pycaster package is currently being hosted on PyPI under this link, while the repo can be found on BitBucket under this link. If you’re interested, I suggest you take a look at the BitBucket repo link first for instructions.
pycaster package really only has a single module, also called
pycaster, and contains a single class called
rayCaster with methods that are simply wrapping the VTK code you saw above. Actually, the VTK code you saw above was pretty much copied from
If you want to jump straight in then installing
pycaster is as easy as:
pip install pycaster
pycaster depends on VTK which, if not already installed, might be a hassle to install through
pip. If you run into trouble with the installation the I suggest you read my first post on Anaconda, follow the instructions to create a nice environment, e.g.,
py27, and install
pycaster as such:
source activate py27 # use 'activate py27' on Windows conda install vtk conda install nose pip install pycaster --no-deps
In this example I will show how to use
pycaster to repeat the above process but this time I’ll be doing so on the STL model of the human skull we extracted in the previous post about surface-extraction. You can find the corresponding notebook here, while the STL model of the skull can be download here. I’ve also repeated the hollow-sphere example we saw above with
pycaster which you can find in this notebook.
The VTK rendering parts of the code we saw before are still the same so I won’t be repeating them here. Please check the notebook for this example.
Loading and rendering the mesh
loadSTL helper-function we used before to… well load the STL, was actually copied from a static method residing under the
rayCaster class with which we initialize a new
rayCaster object. Let’s see how it’s done:
from pycaster import pycaster caster = pycaster.rayCaster.fromSTL("bones.stl", scale=1)
As you can see we initially
pycaster package. Subsequently, we use the
fromSTL static-method of the
rayCaster class to immediately load the STL file, ensure that it was loaded correctly, and then create a new
rayCaster object under the
caster variable. Note that we can also set a
scale for that mesh which is particularly useful when we want to ‘convert’ the mesh units, e.g., from
m etc. If
scale is set to anything but the default
scaleMesh is called during initialization, scales the mesh, and replaces it internally.
If you already have a
vtkPolyDatamesh loaded/created in some other way, you can create a
rayCasterobject through the default constructor as such:
caster = pycaster.rayCaster(mesh)where
vtkPolyDataobject. If you then want to ‘scale’ that mesh as was done above, you can use the
scaleMeshmethod and do so as such:
Using the same VTK code to render the
mesh, which now resides under
caster.mesh, we get the next figure.
You might notice that the skull rendered above displays a much cleaner surface than what we extracted in the ‘Surface Extraction: Creating a mesh from pixel-data using Python and VTK’ post. Well, to be honest I cheated a tad and edited the model through Rhino3D and Geomagic Studio. I did so to place the model around the cartesian origin point, thus eliminating the need to reposition the scene camera, to get rid of all those disjoint meshes, and to clean spikes and holes.
Preparing and rendering the ray
The ‘preparation’ of the ray is exactly the same as in the case of the pure-VTK code. We simply create two
list objects with the coordinates which in this case are:
pSource = [100.0, 100.0, 0.0] pTarget = [0.0, 0.0, 0.0]
We then use the same VTK code as before to render that ray (prior to the actual intersection) resulting in the following figure.
Testing for intersection and extracting the intersection point coordinates
Now here comes why I bothered writing
pycaster. As you saw in the VTK example before, at this point we would need to create and initialize the
vtkOBBTree, create empty
vtkPoints objects, test for intersection, and then convert the coordinates to
tuple objects. However, with
pycaster the whole process is much simpler.
pTarget coordinates we set before, we test for intersection and retrieve the point coordinates using the
castRay method of the
rayCaster class as such:
pointsIntersection = caster.castRay(pSource, pTarget)
which returns a
pointsIntersection. As you can see we saved ourselves a bunch of code and got a clean pythonic object which we can easily
[(62.1171, 62.1171, 0.0000), (56.1171, 56.1171, 0.0000)]
Note that the
castRaymethod will return an empty
listif no intersection points were found and give an
infomessage through the
Using the same VTK rendering approach to display the intersection points as last time we get the following figure.
I should mention another method within the
rayCaster class which IMHO is very useful. The method is called
calcDistanceInSolid and its called in the same way as
castRay, i.e., by providing the ray coordinates.
What this method does is calculate the distance a given ray ‘travels’ within the closed surface, i.e., the solid. While that sounds simple its actually rather complicated as VTK only provides a list of intersection points without specifying whether they’re entry or exit points. Combined with the fact that the ‘source’ and/or ‘target’ of a ray may already reside within the surface, this calculation becomes rather convoluted.
Calling this method in the case of the skull mesh as
caster.calcDistanceInSolid(pSource, pTarget) gives us a distance of ca.
8.485 which you can verify by calculating the euclidian distance between the two intersection points we saw before.
calcDistanceInSolidmethod calculates the total distance ‘travelled’ by a ray and will take into account all intersection points found, i.e., multiple entries and exits. In addition, it will account for the the ‘source’ and ‘target’ points of the ray regardless of whether they were within the surface to begin with or not.
Links & Resources
As my posts are becoming ludicrously long, I decided to start re-listing the links to the various material here so you won’t have to go through the whole post to find them:
Here’s the material used in this post:
- Hollow-sphere example with VTK: IPython Notebook and STL file.
- Hollow-sphere example with
pycaster: IPython Notebook and STL file.
- Human-skull example with
pycaster: IPython Notebook and STL file.
Here are the pages where you can find the
- PyPI page: https://pypi.python.org/pypi/pycaster
- BitBucket repo: https://bitbucket.org/somada141/pycaster
Check these past posts which were used and referenced today:
- Anaconda: The crème de la crème of Python distros
- IPython Notebook & VTK
- Surface Extraction: Creating a mesh from pixel-data using Python and VTK
Don’t forget: all material I’m presenting in this blog can be found under the PyScience BitBucket repository.
Thus concludes another humble post. I hope you enjoyed the whole ray-casting experience as it can be the cornerstone to a lot of super-interesting ray-tracing projects. Off the top of my head you can write code for ray-tracing renderers, ray-tracing physics simulations in optics and acoustics, while ray-casting and ray-tracing is heavily used in game-development (collision detection, rendering, etc).
However, there’s a fair amount of work/code required to get from this example to a ray-tracer as you need to find the mesh-cells where intersection happens, calculate their normal vectors, calculate the reflected/diffracted ray vectors, etc etc. I’m going to show you a lot of this functionality next week in a post dedicated to ray-tracing with VTK.
The material presented here has been tested extensively so you shouldn’t have any trouble reproducing the whole thing. However, if you do then feel free to drop me a comment here and I’ll try to help.
Thanks for reading!