In this post I will show how to use VTK to trace rays emanating from the cell-centers of a source mesh, intersecting with another target mesh, and then show you how to cast subsequent rays bouncing off the target adhering to physics laws. This will include calculating the cell-centers in a mesh, calculating the normal vectors at those cells, vector visualization through glyphs, as well as other elements of visualization like textures and scene-lighting.
In my last post, I used Python and VTK to show you how to perform ray-casting, i.e., intersection tests between arbitrary lines/rays and a mesh, and extraction of the intersection point coordinates through the
Today I’ll take the lessons learned about ray-casting with Python and VTK, and give you the tools to write your own ray-tracing algorithm using Python and VTK. Before proceeding, I strongly recommend that you read the ray-casting post, cause I will be re-using a lot of the functionality presented prior but I won’t repeat myself in too much detail.
We will start by creating the ‘environment’, i.e., the scene, which will comprise a yellowish half-sphere dubbed the
sun, which will act as the ray-source, and a larger nicely textured sphere called
earth which will be the target of those rays.
We will calculate the cell centers of the
sun mesh and cast rays along the directions of the normal vector at each one of those cells. We will then use the
vtkOBBTree functionality we presented in the last post about ray-casting to detect which of these rays intersect with
earth, calculate the appropriate reflected vectors based on the
earth normal vectors, and cast subsequent rays.
Now let me be clear, the code that will be presented today can not, by any stretch of the imagination, be called a fully-fledged ray-tracer. However, I will be presenting all necessary tools you would need to write your own ray-tracing algorithms using Python and VTK.
Ray-Tracing with Python & VTK
You can find the entire IPython Notebook I’ll be presenting today here. It contains a fair bit of code but it was structured in the same way as this post so you can easily look up the different parts of code and see what they do in detail.
As we’ve been doing in the past few posts we need to import
import vtk import numpy
I know I’ve been saying it in nearly every post but here goes again: If for any reason you haven’t managed to get yourself a nice Python installation with
vtk, which are needed for this post, do yourselves a favor and use Anaconda as instructed in this previous post. Alternatively, you can opt for the Enthought Python Distribution (EPD) or Enthought Canopy. Of course there’s a truckload of other distros but I can vouch for the above and they all come with pre-compiled builds of VTK.
The following ‘helper-functions’ are defined at the beginning of today’s notebook and used throughout:
vtk_show(renderer, width=400, height=300): This function allows me to pass a
vtkRendererobject and get a PNG image output of that render, compatible with the IPython Notebook cell output. This code was presented in this past post about VTK integration with an IPython Notebook.
addLine(renderer, p1, p2, color=[0.0, 0.0, 1.0], opacity=1.0): This function uses
vtkto add a line, defined by coordinates under
p2, to a
renderer. In addition, it allows for the color and opacity level of the line to be defined and it was first presented in the previous post about ray-casting.
addPoint(renderer, p, color=[0.0, 0.0, 0.0], radius=0.5): This function uses
vtkto add a point-sized sphere (
0.5) defined by center coordinates under
p. This function was first presented in the previous post about ray-casting, while similar code was detailed in this past post about VTK integration with an IPython Notebook.
l2n = lambda l: numpy.array(l)and
n2l = lambda n: list(n): These are just two simple
lambdafunctions meant to quickly convert a
numpy.ndarrayand vice-versa. The reason I wrote those, is that while VTK returns data like coordinates and vectors in
listobjects, which we need to ‘convert’ to
numpy.ndarrayobjects in order to perform some basic vector math as we’ll see later on. In addition, we often need to feed such vectors and point back to VTK, so back-conversion is often necessary.
As you will see later we’re defining more ‘auxiliary functions’ which haven’t been presented prior to this post so we’re gonna be looking closely into what they do. The above ‘helper-functions’ all contain code that has been presented before at one time or another.
As the code in today’s post deals with a lot of rendering and graphics-related parameters, I decided to set all these parameters as ‘options’ at the beginning of the notebook. You can change any of these and re-run the notebook to ascertain their effect.
You can see these options below but you do not have to pay too much attention to them right now. Apart from a few basic parameters defining attributes of the scene’s objects, they mostly deal with colors. In addition I will be referring back to these options later on while presenting the code.
# SUN OPTIONS # Radius of the sun half-sphere RadiusSun = 10.0 # Distance of sun's center from (0,0,0) DistanceSun = 50.0 # Phi & Theta Resolution of sun ResolutionSun = 6 # EARTH OPTIONS # Radius of the earth sphere RadiusEarth = 150.0 # Phi & Theta Resolution of earth ResolutionEarth = 120 # RAY OPTIONS # Length of rays cast from the sun. Since the rays we # cast are finite lines we set a length appropriate to the scene. # CAUTION: If your rays are too short they won't hit the earth # and ray-tracing won't be possible. Its better for the rays to be # longer than necessary than the other way around RayCastLength = 500.0 # COLOR OPTIONS # Color of the sun half-sphere's surface ColorSun = [1.0, 1.0, 0.0] # Color of the sun half-sphere's edges ColorSunEdge = [0.0, 0.0, 0.0] # Color of the earth sphere's edges ColorEarthEdge = [1.0, 1.0, 1.0] # Background color of the scene ColorBackground = [0.0, 0.0, 0.0] # Ambient color light ColorLight = [1.0, 1.0, 0.0] # Color of the sun's cell-center points ColorSunPoints = [1.0, 1.0, 0.0] # Color of the sun's cell-center normal-vector glyphs ColorSunGlyphs = [1.0, 1.0, 0.0] # Color of sun rays that intersect with earth ColorRayHit = [1.0, 1.0, 0.0] # Color of sun rays that miss the earth ColorRayMiss = [1.0, 1.0, 1.0] # Opacity of sun rays that miss the earth OpacityRayMiss = 0.5 # Color of rays showing the normals from points on earth hit by sun rays ColorEarthGlyphs = [0.0, 0.0, 1.0] # Color of sun rays bouncing off earth ColorRayReflected = [1.0, 1.0, 0.0]
We will start by creating the ‘environment’ within which we’ll perform our ray-tracing. This environment will comprise a half-sphere dubbed
sun from which we will cast rays. Another sphere, which we will call
earth, will receive and reflect those rays at appropriate angles.
Note that the
earthobjects are by no means ‘in-scale’. Far from it actually (whole thing kinda looks like an earth-ball with a ceiling light above). I just named them as such to provide an apt analogy to what we’re doing here.
We start by creating the
sun half-sphere through
# Create and configure then sun half-sphere sun = vtk.vtkSphereSource() sun.SetCenter(0.0, DistanceSun, 0.0) sun.SetRadius(RadiusSun) sun.SetThetaResolution(ResolutionSun) sun.SetPhiResolution(ResolutionSun) sun.SetStartTheta(180) # create a half-sphere
source->mapper->actor process to create a sphere in VTK has been shown time after time and you can read about the mechanics in this previous post. However, I’ll go over a few novelties here. As you can see, upon creating a new
sun we set several parameters. The center coordinates and radius are defined through the
RadiusSun options discussed in the Options section.
What’s really important to note here is the
phi resolutions of the
sun sphere set through the
SetPhiResolution methods respectively. VTK creates these spheres through the use of spherical coordinates and these resolution parameters define the number of points along the longitude and latitude directions respectively. A lower resolution will result in less points and therefore less triangles defining the sphere, making it look like something of a rough polygon.
However, as I mentioned in the Summary we will be casting a ray per triangle of the
sun mesh. Hence, in order to keep the scene ‘clean’ and lower the computational cost, we set this resolutions to a mere
6 in the Options. Nonetheless, feel free to change this value in order to cast as many or as few rays as you wish. I’ll be showing the result of a resolution of
20 at the end of this post.
A last point I’d like to make is the usage of the
SetStartTheta method. As I just said, VTK uses spherical coordinates to ‘design’ these spherical objects. Through the
SetStopPhi methods we define the starting and stopping angles (in degrees) for this sphere, allowing us to create spherical segments instead of a full sphere. In this example, and in the interest of keeping the scene ‘lean’ we’re only creating a ‘half-sphere’ through
# Create mapper mapperSun = vtk.vtkPolyDataMapper() mapperSun.SetInput(sun.GetOutput()) # Create actor actorSun = vtk.vtkActor() actorSun.SetMapper(mapperSun) actorSun.GetProperty().SetColor(ColorSun) #set color to yellow actorSun.GetProperty().EdgeVisibilityOn() # show edges/wireframe actorSun.GetProperty().SetEdgeColor(ColorSunEdge) #render edges as white
Once more, check this previous post to see what the deal with mappers and actors is in VTK. What I should mention here is that apart from using the
GetProperty() method to access the properties of the
actorSun object and set its color to
ColorSun (set under Options), we also do so to enable the visibility of that mesh’s edges and set them to a color of
ColorSunEdge. This way we’re visualizing the ‘wireframe’ of that object in order to see the individual mesh cells which come into play later.
Finally as we’ve done in all posts dealing with VTK we create a new
vtkRenderer, add the
actorSun object, set some camera properties, and use the
vtk_show ‘helper-function’ to render the scene:
renderer = vtk.vtkRenderer() renderer.AddActor(actorSun) renderer.SetBackground(ColorBackground) # Modify the camera with properties defined manually in ParaView camera = renderer.MakeCamera() camera.SetPosition(RadiusEarth, DistanceSun, RadiusEarth) camera.SetFocalPoint(0.0, 0.0, 0.0) camera.SetViewAngle(30.0) renderer.SetActiveCamera(camera) vtk_show(renderer, 600, 600)
Take a look at this past post about VTK integration with an IPython Notebook for an explanation of the
vtkRenderer class, and the post about surface extraction to see what gives with the camera. What’s important to note is that we’ve created a
vtkRenderer object under
renderer to which we will continue adding new
vtkActor objects as we go, thus enriching the scene.
The result of
vtk_show can be seen in the following figure:
Here you can see the effect of that ‘resolution’ stuff I was talking about before. The
sunhalf-sphere is jagged and rough. Nonetheless, it still comprises plenty of cells and therefore ray-sources for our example.
Let’s go onto creating
earth. As during my posts I’ve showed you how to create spheres a trillion times I thought I’d kick it up a notch and show you how to texture one, giving it a bit of razzle-dazzle :). We start by creating a new sphere under
earth as such:
# Create and configure the earth sphere earth = vtk.vtkSphereSource() earth.SetCenter(0.0, -RadiusEarth, 0.0) earth.SetThetaResolution(ResolutionEarth) earth.SetPhiResolution(ResolutionEarth) earth.SetRadius(RadiusEarth)
Typical stuff as you can see, just keep in mind that the
earth variable holds the pointer to the, now configured,
vtkSphereSource object. You’ll also note that all parameters are set to values defined in the Options section so take a look if you will.
Now let’s move on to texturing. I used a texture of the earth I downloaded as a JPEG file from http://planetpixelemporium.com/download/download.php?earthmap1k.jpg and placed along this post’s notebook. Alternatively, you can download it from the blog repository under here. Firstly, we need to load that image as such:
# Load a JPEG file with an 'earth' texture downloaded from the above link reader = vtk.vtkJPEGReader() reader.SetFileName("earthmap1k.jpg") reader.Update()
All we do is use the
vtkJPEGReader class to create
reader, set the filename, and read in the image through
I want to specify, that as the
readeris not used directly but is rather a part of the VTK pipeline, we don’t really need to call
Update. That method will be called by subsequent objects that use
readeras an input. However, while debugging VTK code, progressively updating the pipeline allows us to verify which part of the pipeline didn’t work instead of getting an arbitrary error at some point down the road.
Next, we create the actual texture:
# Create a new 'vtkTexture' and set the loaded JPEG texture = vtk.vtkTexture() texture.SetInputConnection(reader.GetOutputPort()) texture.Update()
Here we create a new
vtkTexture object and, using the
GetOutputPort methods, connect its input to the output of the
reader object holding the texture image.
So far so good right? What needs to be done now to complete the texturing, is to map this texture to our
earth sphere. Now this might be a tad confusing so I’ll break it down:
# Map the earth texture to the earth sphere map_to_sphere = vtk.vtkTextureMapToSphere() map_to_sphere.SetInputConnection(earth.GetOutputPort()) map_to_sphere.PreventSeamOn() texture.Update()
As a first step we create a new
vtkTextureMapToSphere object under
map_to_sphere which “generates 2D texture coordinates by mapping input dataset points onto a sphere”. Check the code carefully, and you’ll realize we’re not connecting
texture but only with
earth! This class will just calculate the coordinates based on the
earth sphere but does not set a texture yet (that occurs later).
And now for the finishing touches:
# Create a new mapper with the mapped texture and sphere mapperEarth = vtk.vtkPolyDataMapper() mapperEarth.SetInputConnection(map_to_sphere.GetOutputPort()) # Create actor actorEarth = vtk.vtkActor() actorEarth.SetMapper(mapperEarth) actorEarth.SetTexture(texture) actorEarth.GetProperty().EdgeVisibilityOn() # show edges/wireframe actorEarth.GetProperty().SetEdgeColor(ColorEarthEdge) #render edges as white renderer.AddActor(actorEarth) vtk_show(renderer, 600, 600)
You should be familiar with the
vtkActor classes so I won’t repeat myself. What you should pay attention to is that instead of using
earth as a source for
mapperEarth we instead use
map_to_sphere which contains the
vtkTextureMapToSphere object we just saw! In addition, upon creating
actorEarth we finally set the
actorEarth.SetTexture(texture). As you can see the
texture goes straight into the appropriate
vtkActor object, which along with the texture-mapping provides us with a nicely textured sphere.
Finally, we make the
earth wireframe visible with a
ColorEarthEdge color, add
renderer, and use
vtk_show to render the scene resulting in the following figure:
Adding lighting to the scene
Now you might say “what gives? we went to so much trouble to texture that stupid ball and it looks all dark and crummy!”. Well at least I obviously thought so, thus I decided to shed a little light on the situation (I know, stupidest pun ever).
Before I continue, I want to mention that all VTK scenes comes with a default ‘headlight’ which follows the camera and illuminates the scene. In most cases that’s enough to see our rendering but often we need a lil’ extra.
In this case I added a light point-source through the
vtkLight class. Let’s see how that was done:
# Create a new vtkLight light = vtk.vtkLight() # Set its ambient color to yellow light.SetAmbientColor(ColorLight) # Set a 180 degree cone angle to avoid a spotlight-effect light.SetConeAngle(180) # Set its position to the sun's center light.SetPosition(sun.GetCenter()) # Set its focal-point to the earth's center light.SetFocalPoint(earth.GetCenter()) # Set it as part of the scene (positional) and not following the camera light.SetPositional(True) renderer.AddLight(light) vtk_show(renderer, 600, 600)
Initially, we create a
vtkLight object under
light, and set an appropriate
ColorLight which was defined at the Options section. As you can see I set the position of the light-source to the center of the
sun half-sphere, which just makes sense (you know,
sun and all that), through
light.SetPosition(sun.GetCenter()), while its focal point is the center of the
What’s important to note here is that these lights are typically ‘spotlights’ adding a ‘narrow’ beam of light to the scene. By using
SetConeAngle(180) we create a uniform semi-spherical light that evenly illuminates the entire scene residing beneath the center of the
sun half-sphere. Lastly, please pay attention to the
SetPositional(True) method which makes this light a constant part of the scene, forcing it in place, and not allowing it to move with respect to the scene’s camera.
light to the
renderer, we once more use
vtk_show to render the scene and get this figure (makes quite the difference right?):
‘Prepare’ the Sun Rays
As I said in the Summary, we will be casting rays from each cell-center of the
sun mesh following the direction of the normal vectors of those cells. Naturally, we first need to calculate these quantities, thus ‘defining’ the rays.
Calculate the cell-centers of the sun half-sphere
Firstly, we need to calculate the coordinates at the center of each cell on the
sun mesh. Thankfully VTK provides us with a nifty class called
vtkCellCenters to do so:
cellCenterCalcSun = vtk.vtkCellCenters() cellCenterCalcSun.SetInput(sun.GetOutput()) cellCenterCalcSun.Update()
As you can see this is as simple as VKT gets. We merely create a
vtkCellCenters object under
cellCenterCalcSun, and connect its input with
sun.GetOutput(). Upon using
Update, the calculation is complete.
Now let us visualize those points using the
vtk_show helper functions:
# Get the point centers from 'cellCenterCalc' pointsCellCentersSun = cellCenterCalcSun.GetOutput(0) # Loop through all point centers and add a point-actor through 'addPoint' for idx in range(pointsCellCentersSun.GetNumberOfPoints()): addPoint(renderer, pointsCellCentersSun.GetPoint(idx), ColorSunPoints) vtk_show(renderer, 600, 600)
We first ‘extract’ the cell-centers through
cellCenterCalcSun.GetOutput(0), storing the result of
vtkPoints type under
pointsCellCentersSun. Subsequently, we loop through these points and use
addPoint to add them to the
vtkRenderer object we created before and which now resides under
renderer. Note that we’re looping using the
range built-in and the
GetNumberOfPoints() method to get the total number of points found. The resulting figure can then be seen below.
Calculate normal vectors at the center of each cell
Now that we have the cell-centers, which will act as the ‘source points’ of the
sun rays, we need to calculate the normal vectors at those points which will define the directions of the rays. Once more, VTK provides the tools but doesn’t make it easy or clear for us (you wouldn’t appreciate it working if it was easy, would you now 🙂 ?). Let’s see how it’s done:
# Create a new 'vtkPolyDataNormals' and connect to the 'sun' half-sphere normalsCalcSun = vtk.vtkPolyDataNormals() normalsCalcSun.SetInputConnection(sun.GetOutputPort()) # Disable normal calculation at cell vertices normalsCalcSun.ComputePointNormalsOff() # Enable normal calculation at cell centers normalsCalcSun.ComputeCellNormalsOn() # Disable splitting of sharp edges normalsCalcSun.SplittingOff() # Disable global flipping of normal orientation normalsCalcSun.FlipNormalsOff() # Enable automatic determination of correct normal orientation normalsCalcSun.AutoOrientNormalsOn() # Perform calculation normalsCalcSun.Update()
So we start by creating a new
vtkPolyDataNormals object under
normalsCalcSun (remember this variable name, we’ll be using it lots later on). We then ‘connect’ this object to the
sun mesh where we want those normals to be calculated.
The rest boils down to configuring the
vtkPolyDataNormals class to give us what we want. You can see what each call does in the comments above but I should stress a few points here.
vtkPolyDataNormals can calculate the normal vectors at the mesh points and/or cells. However, we only want the latter, so we turn off calculation at points with
ComputePointNormalsOff() and only enable calculation at cells with
Subsequently, we turn ‘splitting’ off through
SplittingOff(). First of all that only makes sense when calculating point-normals. What it would do is ‘split’, thus create, multiple normals at points belonging to cells with very sharp edges, i.e., steep angles between cells. However, we only want cell normals so we don’t care about it too much.
We then want to make sure that the normals have a correct orientation, i.e., that they would be ‘pointing’ outwards and not towards the center of the
sun half-sphere. To ensure that, we first disable global flipping through ‘FlipNormalsOff()’, and we enable automatic orientation determination through
AutoOrientNormalsOn(). Now this last call is a god-send as it will make sure the normals point outwards which we sorely need to correctly cast our rays.
However, I want to point you to the
vtkPolyDataNormals docs, and in particular the docstring for the
AutoOrientNormalsOn method which reads: “Turn on/off the automatic determination of correct normal orientation. NOTE: This assumes a completely closed surface (i.e. no boundary edges) and no non-manifold edges. If these constraints do not hold, all bets are off. This option adds some computational complexity, and is useful if you don’t want to have to inspect the rendered image to determine whether to turn on the FlipNormals flag. However, this flag can work with the FlipNormals flag, and if both are set, all the normals in the output will point ”inward“.”
What that means is that while
AutoOrientNormalsOn is crazy-useful to ensure correct orientation, it comes with stringent requirements, and its not always guaranteed to succeed. Thankfully, our
sun mesh fits these criteria and calling
Update() completes the calculation.
Visualize normal vectors at the cell-centers of sun’s surface as glyphs
Before proceeding on to casting and tracing the
sun rays we will visualize the normal vectors calculated on the
sun mesh as glyphs using the
Before I show you the code lets take a look at the docstring of the
vtkGlyph3D class: “vtkGlyph3D is a filter that copies a geometric representation (called a glyph) to every point in the input dataset. The glyph is defined with polygonal data from a source filter input. The glyph may be oriented along the input vectors or normals, and it may be scaled according to scalar data or vector magnitude”.
What we’re going to do here, is create a single arrow through the
vtkArrowSource class and use it as a glyph. Then we’ll use the normal vectors we calculated before, and which are stored under
normalsCalcSun, to place and orient those glyphs. Let’s inspect the code:
# Create a 'dummy' 'vtkCellCenters' to force the glyphs to the cell-centers dummy_cellCenterCalcSun = vtk.vtkCellCenters() dummy_cellCenterCalcSun.VertexCellsOn() dummy_cellCenterCalcSun.SetInputConnection(normalsCalcSun.GetOutputPort()) # Create a new 'default' arrow to use as a glyph arrow = vtk.vtkArrowSource() # Create a new 'vtkGlyph3D' glyphSun = vtk.vtkGlyph3D() # Set its 'input' as the cell-center normals calculated at the sun's cells glyphSun.SetInputConnection(dummy_cellCenterCalcSun.GetOutputPort()) # Set its 'source', i.e., the glyph object, as the 'arrow' glyphSun.SetSourceConnection(arrow.GetOutputPort()) # Enforce usage of normals for orientation glyphSun.SetVectorModeToUseNormal() # Set scale for the arrow object glyphSun.SetScaleFactor(5) # Create a mapper for all the arrow-glyphs glyphMapperSun = vtk.vtkPolyDataMapper() glyphMapperSun.SetInputConnection(glyphSun.GetOutputPort()) # Create an actor for the arrow-glyphs glyphActorSun = vtk.vtkActor() glyphActorSun.SetMapper(glyphMapperSun) glyphActorSun.GetProperty().SetColor(ColorSunGlyphs) # Add actor renderer.AddActor(glyphActorSun) vtk_show(renderer, 600, 600)
Now this code is a little convoluted so I’ll break it down. The first very important thing to mention is that we do not simply use the normals calculated beforehand stored within
normalsCalcSun. As you can see in the first lines of the snippet above, we feed those normals to a ‘dummy’
dummy_cellCenterCalcSun. This forces the glyphs to be placed on the cell-centers of the
Subsequently, we create the glyphs. Note that we first create
vtkArrowSource object representing a default arrow, which we will use as the ‘base’ glyph. Now this glyph can be any ‘source’ class, e.g. a cone through the
vtkConeSource class, a sphere through the
vtkSphereSource class, etc etc. The only reason I chose an arrow was cause it nicely shows the direction the
sun rays will follow.
We then create a new
vtkGlyph3D object under
glyphSun. The important thing to note here is the difference between the
SetSourceConnection methods. The latter just connects to the ‘source’ glyph, i.e., the
arrow in our case. The
SetInputConnection call is given
dummy_cellCenterCalcSun.GetOutputPort(), i.e., the normal vectors calculated and then positioned at the cell-centers of the
Henceforth things are simple: we enforce orientation of the created glyphs to the supplied normal vectors through
SetVectorModeToUseNormal(), and we provide a ‘scale factor’ for the base glyphs, which will just uniformly scale the
arrow and its default size. The rest has been shown a thousand times: we create a
vtkPolyDataMapper to map the create object to graphics primitives and connect to the
glyphSun output. We then create a standard
vtkActor, connect to the aforementioned mapper, and use its
GetProperty() method to
Finally, using the
vtk_show helper-function yields the following figure. As you can see we have visualized all normal vectors with arrows, showing the direction our
sun rays will follow.
Prepare for ray-tracing
We’re finally getting to the ray-tracing part of the post. All we now need to do is prepare the
vtkOBBTree object for
earth as I showed in the last post on ray-casting. If you haven’t read it then I strongly recommend that you do now cause I won’t explain the details again.
Firstly, we create a new
vtkOBBTree object with the
earth mesh where we’re going to test for intersection with rays coming from the
sun. This is done as simply as this:
obbEarth = vtk.vtkOBBTree() obbEarth.SetDataSet(earth.GetOutput()) obbEarth.BuildLocator()
obbEarth in mind cause we’ll be using extensively to test for intersections later.
Then we need calculate the normal vectors at all cell-centers on the
earth mesh as we’re going to need that information to cast reflected rays. The process here is exactly the same as the one we used to calculate the normal vectors for the
# Create a new 'vtkPolyDataNormals' and connect to the 'earth' sphere normalsCalcEarth = vtk.vtkPolyDataNormals() normalsCalcEarth.SetInputConnection(earth.GetOutputPort()) # Disable normal calculation at cell vertices normalsCalcEarth.ComputePointNormalsOff() # Enable normal calculation at cell centers normalsCalcEarth.ComputeCellNormalsOn() # Disable splitting of sharp edges normalsCalcEarth.SplittingOff() # Disable global flipping of normal orientation normalsCalcEarth.FlipNormalsOff() # Enable automatic determination of correct normal orientation normalsCalcEarth.AutoOrientNormalsOn() # Perform calculation normalsCalcEarth.Update()
Just remember that
normalsCalcEarth now holds the normal vectors for the
Finally, we define two ‘auxiliary-functions’ to nicely wrap the intersection testing functionality offered by the
def isHit(obbTree, pSource, pTarget): r"""Returns True if the line intersects with the mesh in 'obbTree'""" code = obbTree.IntersectWithLine(pSource, pTarget, None, None) if code==0: return False return True def GetIntersect(obbTree, pSource, pTarget): # Create an empty 'vtkPoints' object to store the intersection point coordinates points = vtk.vtkPoints() # Create an empty 'vtkIdList' object to store the ids of the cells that intersect # with the cast rays cellIds = vtk.vtkIdList() # Perform intersection code = obbTree.IntersectWithLine(pSource, pTarget, points, cellIds) # Get point-data pointData = points.GetData() # Get number of intersection points found noPoints = pointData.GetNumberOfTuples() # Get number of intersected cell ids noIds = cellIds.GetNumberOfIds() assert (noPoints == noIds) # Loop through the found points and cells and store # them in lists pointsInter =  cellIdsInter =  for idx in range(noPoints): pointsInter.append(pointData.GetTuple3(idx)) cellIdsInter.append(cellIds.GetId(idx)) return pointsInter, cellIdsInter
Again, if you haven’t read the last post on ray-casting, do so and the above will make complete sense to you.
isHit function will simply return
False depending on whether a given ray intersects with
obbTree, which in our case will only be
The second function is a little more complex but the mechanics were detailed in the previous post. I’ve included comments in the code so you can get what it does but I’m not going to spoon-feed this to you again. Read up :).
In a nutshell it will return two
cellIdsInter. The former will contain a series of
tuple objects with the coordinates of the intersection points. The latter will contain the ‘id’ of the mesh cells that were ‘hit’ by that ray. This information is vital as through these ids we’ll be able to get the correct normal vector for that
earth cell and calculate the appropriate reflected vector as we’ll see below.
Perform ray-casting operations and visualize different aspects
Finally the moment you’ve all been waiting for! The ray-tracing! Now since the code to perform the whole operation is too much to explain in a single go, I decided to repeat the process three times, visualizing and explaining different aspects of it every time.
Cast rays, test for intersections, and visualize the rays that hit
In this first step we will only cast the rays from the
sun, test for their intersection, or lack thereof, and render the rays that do hit the
earth with a different color than the ones that miss it. In addition, we’ll also render points at those intersection points as they will be the ‘source’ points for the reflected rays later on.
The code is the following:
# Extract the normal-vector data at the sun's cells normalsSun = normalsCalcSun.GetOutput().GetCellData().GetNormals() # Loop through all of sun's cell-centers for idx in range(pointsCellCentersSun.GetNumberOfPoints()): # Get coordinates of sun's cell center pointSun = pointsCellCentersSun.GetPoint(idx) # Get normal vector at that cell normalSun = normalsSun.GetTuple(idx) # Calculate the 'target' of the ray based on 'RayCastLength' pointRayTarget = n2l(l2n(pointSun) + RayCastLength*l2n(normalSun)) # Check if there are any intersections for the given ray if isHit(obbEarth, pointSun, pointRayTarget): # intersections were found # Retrieve coordinates of intersection points and intersected cell ids pointsInter, cellIdsInter = GetIntersect(obbEarth, pointSun, pointRayTarget) # Render lines/rays emanating from the sun. Rays that intersect are # rendered with a 'ColorRayHit' color addLine(renderer, pointSun, pointsInter, ColorRayHit) # Render intersection points addPoint(renderer, pointsInter, ColorRayHit) else: # Rays that miss the earth are rendered with a 'ColorRayMiss' colors # and a 25% opacity addLine(renderer, pointSun, pointRayTarget, ColorRayMiss, opacity=0.25) vtk_show(renderer, 600, 600)
As you no doubt remember, cause you’ve been paying attention all this time, the
normalsCalcSun is a
vtkPolyDataNormals object which holds the normal vectors calculated at the cell-centers of the
sun mesh. As you can see from the 1st line of the code, retrieving the actual normal vector data is easy but not intuitive. First we get all output out through
GetOutput. However, the output could have point-data and/or cell-data depending on how we configured the class. In our case we want the cell-data which we retrieve through
GetCellData() followed by
GetNormals() which will give a set of normal vectors under
Afterwards, as we want to cast a ray from every cell-center on the
sun mesh, we loop through these points stored under
idx variable will hold an index to that point.
We store the coordinates of each such point under
pointsCellCentersSun.GetPoint(idx), while we store the corresponding normal vector under
normalsSun.GetTuple(idx). Now we have the ‘source’ point and direction our ray should follow.
In our example, however, rays are merely lines of finite length. Due to that we have to define a ‘target’ point which we want to make sure it won’t miss the
earth due to an insufficient length. To that end I’ve defined a variable
RayCastLength under the Options, which naturally will define the length of the ‘ray’. As such we use the following line (repeating it here) to calculate the coordinates of the ‘target’ point:
pointRayTarget = n2l(l2n(pointSun) + RayCastLength*l2n(normalSun))
Note that we use the
l2n helper-functions to quickly convert a
tuple object as retrieved from VTK, into a
numpy.ndarray object and vice-versa. As I mentioned in the Helper-Functions section, this conversion allows us to use
numpy and perform some basic vector-math, which in this case is vector addition and multiplication.
Now we finally have all we need to cast that pesky ray, i.e. the ‘source’ and ‘target’ coordinates! We first use the
isHit auxiliary-function to cast that ray and see if it intersects with
earth as such:
isHit(obbEarth, pointSun, pointRayTarget)
If we get a
False back then the ray does not intersect and its not worth processing further but we render it anyway through
addLine with a
ColorRayMiss color, and a 25% opacity (in our case resulting in a faded whitish ray).
However, should the ray intersect with
earth then its worth processing. We then re-cast that ray and retrieve the intersection point coordinates and ids of the
earth cells the ray hit through the
GetIntersect auxiliary-function we defined before (repeating the line here):
pointsInter, cellIdsInter = GetIntersect(obbEarth, pointSun, pointRayTarget)
As I explained before, the
pointsInter variable now holds a
tuple objects with the coordinates of each intersection point in order. The
cellIdsInter contains simple
int entries with the id of the mesh cells that were this but we won’t use that information in this step. What we will however use is
pointsInter, i.e., the 1st point that was hit by the ray.
We consider the
earthto be a reflective body so we don’t care about the following intersection points as the rays that impinge on
earthare supposed to be reflected entirely.
The rest is super-simple. We just use
addPoint to render the ray that hits the earth and the 1st intersection point, which under default Options should appear as yellow.
After looping through all cell-centers in the
sun, cast those rays, and define whether they hit or miss, we finally call
vtk_show and end up with the next figure.
Visualize the normal vectors at the cells where sun’s rays intersect with earth
During this next step we’ll add a lil’ more complexity, and visualize the
earth normal vectors at the points where rays from the
sun intersect and we’ll do so again through glyphs. However, unlike the case where we rendered all normals on the
sun surface, here we only want to render the normals at points were intersections were detected. Thus, the glyph process is slightly different and allows for more flexibility.
# Extract the normal-vector data at the sun's cells normalsSun = normalsCalcSun.GetOutput().GetCellData().GetNormals() # Extract the normal-vector data at the earth's cells normalsEarth = normalsCalcEarth.GetOutput().GetCellData().GetNormals() # Create a dummy 'vtkPoints' to act as a container for the point coordinates # where intersections are found dummy_points = vtk.vtkPoints() # Create a dummy 'vtkDoubleArray' to act as a container for the normal # vectors where intersections are found dummy_vectors = vtk.vtkDoubleArray() dummy_vectors.SetNumberOfComponents(3) # Create a dummy 'vtkPolyData' to store points and normals dummy_polydata = vtk.vtkPolyData() # Loop through all of sun's cell-centers for idx in range(pointsCellCentersSun.GetNumberOfPoints()): # Get coordinates of sun's cell center pointSun = pointsCellCentersSun.GetPoint(idx) # Get normal vector at that cell normalSun = normalsSun.GetTuple(idx) # Calculate the 'target' of the ray based on 'RayCastLength' pointRayTarget = n2l(l2n(pointSun) + RayCastLength*l2n(normalSun)) # Check if there are any intersections for the given ray if isHit(obbEarth, pointSun, pointRayTarget): # Retrieve coordinates of intersection points and intersected cell ids pointsInter, cellIdsInter = GetIntersect(obbEarth, pointSun, pointRayTarget) # Get the normal vector at the earth cell that intersected with the ray normalEarth = normalsEarth.GetTuple(cellIdsInter) # Insert the coordinates of the intersection point in the dummy container dummy_points.InsertNextPoint(pointsInter) # Insert the normal vector of the intersection cell in the dummy container dummy_vectors.InsertNextTuple(normalEarth) # Assign the dummy points to the dummy polydata dummy_polydata.SetPoints(dummy_points) # Assign the dummy vectors to the dummy polydata dummy_polydata.GetPointData().SetNormals(dummy_vectors) # Visualize normals as done previously but using # the 'dummyPolyData' arrow = vtk.vtkArrowSource() glyphEarth = vtk.vtkGlyph3D() glyphEarth.SetInput(dummy_polydata) glyphEarth.SetSourceConnection(arrow.GetOutputPort()) glyphEarth.SetVectorModeToUseNormal() glyphEarth.SetScaleFactor(5) glyphMapperEarth = vtk.vtkPolyDataMapper() glyphMapperEarth.SetInputConnection(glyphEarth.GetOutputPort()) glyphActorEarth = vtk.vtkActor() glyphActorEarth.SetMapper(glyphMapperEarth) glyphActorEarth.GetProperty().SetColor(ColorEarthGlyphs) renderer.AddActor(glyphActorEarth) vtk_show(renderer, 600, 600)
Firstly, we extract the normal vectors for the
earth, much as we did before, and store them under
Then, we create two ‘dummy’ containers that will store the intersection point coordinates (
dummy_points) and the
earth normals (
sun rays intersected. We will use this information to render those, and only those, normals on
earth. Note that
dummy_points is of type
dummy_vectors is of
vtkDoubleArray type. We ‘configure’
dummy_vectors to have
3 components through
SetNumberOfComponents which are going to be the vector’s
z, components. Finally, we define
dummy_polydata of type
vtkPolyData, which will ‘house’ the ‘dummy’ points and vectors and which we will pass to the
vtkGlyph3D class later on.
The loop in this step is very similar to the one we presented in the previous step. We simply go through all cell-centers on the
sun mesh, cast-rays, and test for their intersection. The key difference is that instead of rendering rays, we now use the
cellIdsInter to get the cell id of the 1st cell, i.e., the cell the
sun ray first intersects with. We then use that id to retrieve the normal vector on the
earth where the ray hit as such:
normalEarth = normalsEarth.GetTuple(cellIdsInter)
Subsequently, we use the
InsertNextTuple methods to ‘push’ the retrieved point coordinates and normal vectors into the dummy containers. Once the loop is complete, we set those points and vectors to the
dummy_polydata container as points and normals. We do this through the
GetPointData().SetNormals() methods respectively.
The rest is mostly the same as the previous glyph rendering we did for the
sun normals. The only difference is that the source of the
glyphEarth object is now the
dummy_polydata object we just composed. As you can understand, this approach allowed us to render fully-customized glyphs. The result of this step can be seen in the next figure where we can now see the
earth normal vectors where
sun rays intersect.
Calculate and visualize reflected rays
Here comes the final step. We now have all information we need to cast rays from the
sun, detect which ones hit
earth, and use vector math to cast subsequent rays that are reflected off the
earth surface with the appropriate orientation.
We first define a little ‘auxiliary-function’ called
calcVecR that calculates a reflected vector from the incident and normal vectors:
def calcVecR(vecInc, vecNor): vecInc = l2n(vecInc) vecNor = l2n(vecNor) vecRef = vecInc - 2*numpy.dot(vecInc, vecNor)*vecNor return n2l(vecRef)
As you can see we use the
n2l helper-functions to quickly convert vectors from
tuple objects to
numpy.ndarray objects and vice-versa. A nice article detailing this type of vector math used in ray-tracing can be seen here.
Finally, we repeat the whole process we saw before but this time we only calculate and render the reflected rays:
# Extract the normal-vector data at the sun's cells normalsSun = normalsCalcSun.GetOutput().GetCellData().GetNormals() # Extract the normal-vector data at the earth's cells normalsEarth = normalsCalcEarth.GetOutput().GetCellData().GetNormals() # Loop through all of sun's cell-centers for idx in range(pointsCellCentersSun.GetNumberOfPoints()): # Get coordinates of sun's cell center pointSun = pointsCellCentersSun.GetPoint(idx) # Get normal vector at that cell normalSun = normalsSun.GetTuple(idx) # Calculate the 'target' of the ray based on 'RayCastLength' pointRayTarget = n2l(l2n(pointSun) + RayCastLength*l2n(normalSun)) # Check if there are any intersections for the given ray if isHit(obbEarth, pointSun, pointRayTarget): # Retrieve coordinates of intersection points and intersected cell ids pointsInter, cellIdsInter = GetIntersect(obbEarth, pointSun, pointRayTarget) # Get the normal vector at the earth cell that intersected with the ray normalEarth = normalsEarth.GetTuple(cellIdsInter) # Calculate the incident ray vector vecInc = n2l(l2n(pointRayTarget) - l2n(pointSun)) # Calculate the reflected ray vector vecRef = calcVecR(vecInc, normalEarth) # Calculate the 'target' of the reflected ray based on 'RayCastLength' pointRayReflectedTarget = n2l(l2n(pointsInter) + RayCastLength*l2n(vecRef)) # Render lines/rays bouncing off earth with a 'ColorRayReflected' color addLine(renderer, pointsInter, pointRayReflectedTarget, ColorRayReflected) vtk_show(renderer, 600, 600)
The only ‘new’ lines here to which you should pay attention, reside within the loop and are the following:
# Calculate the incident ray vector vecInc = n2l(l2n(pointRayTarget) - l2n(pointSun)) # Calculate the reflected ray vector vecRef = calcVecR(vecInc, normalEarth) # Calculate the 'target' of the reflected ray based on 'RayCastLength' pointRayReflectedTarget = n2l(l2n(pointsInter) + RayCastLength*l2n(vecRef))
As you can see, we just calculate the vector of the incident ray, a ray which we know hits the
earth, by subtracting the coordinates of the ray’s ‘target’ point on
earth from the ray’s ‘source’ point on the
Then we use the simple auxiliary-function
calcVecR to calculate the vector of the reflected ray through
vecRef = calcVecR(vecInc, normalEarth). Finally we calculate a ‘target’ point for this reflected ray, much as we did when we first cast the rays from the
sun, as such:
pointRayReflectedTarget = n2l(l2n(pointsInter) + RayCastLength*l2n(vecRef))
At long last, we just render these reflected rays through
vtk_show and get the next figure.
As you can see, we just went through all the necessary components to perform ray-tracing with a given ‘source’ object being the
sun, and a ‘target’ object being the
earth. For a full ray-tracer, one would need to perform this sort of calculations and tests for each ‘target’ object and keep casting rays. In a more realistic ray-tracer two additional concepts would have to be accounted for:
- Each ray would need to carry a certain amount of ‘energy’ which would gradually deplete, eventually running out, thus providing a means of terminating an otherwise endless loop.
- The ‘target’ object would not merely reflect the impinging rays but part of the ray’s energy would be reflected while the rest would refract through the ‘target’ object. Therefore, every ray-hit would actually result in two subsequent rays that would need to be cast and traced. However, that would require our different objects to exhibit ‘material’ characteristics, thus defining the reflective and refractive indices as well as attenuation that would further deplete the ray’s energy as it passes through such an object.
However, the above would result in code that too much for a conceptual post and I’ll leave it up to you to extend it for yourselves. I just hope I’ve provided you with all you needed to know to realize something of the sort.
One last thing before I close. I mentioned in the beginning of this post, such a long time ago, that we needed to have a low ‘resolution’ for the
sun mesh in order to get a small number of triangles, and therefore rays. If we had a very refined mesh on the
sun not only would we significantly increase the number of rays, and therefore boggle the graphics, but we’d also be faced with quite a bit of computational weight. Regardless, here’s the final render showing the entire scene if we were to set the
ResolutionSun variable under Options from its default value of
6 to a value of
Pretty gorgeous albeit messy, wouldn’t you say 🙂 ?
Links & Resources
Here’s the material used in this post:
- IPython Notebook showing the entire process.
- Earth Texture.
- Reflections and Refractions in Ray-Tracing by Bram de Greve: Article explaining the math and physics used in ray-tracing.
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
- Ray Casting with Python and VTK: Intersecting lines/rays with surface meshes
Don’t forget: all material I’m presenting in this blog can be found under the PyScience BitBucket repository.
I hope you enjoyed going through this post as much as I enjoyed putting it together. This was by far my longest post yet, and it took a preposterous amount of my free time to prepare, but I hope I managed to introduce you to a large variety of different VTK classes which you can utilize in a gazillion of different ways.
I’ve tested the accompanying notebook on Mac OS X 10.9.5, Windows 8.1, and Linux Mint 17, all running the latest Anaconda Python with a
py27 environment as discussed in my first post. Thus, I have no reason to believe that you’ll run into trouble but if you do then you know the drill (comment here and I’ll get back to you) :).
One last quick note: I know I never mentioned it before, cause – to my dismay – I tend to consider a few things as common knowledge, but VTK has not yet been wrapped for Python 3 even though the porting has been in the works for some time. For better or for worse, Python+VTK code will need to be executed on Python 2 for the time being :).