Figuring Out Stereo Volume Rendering

Last week I was approached by a friend and collaborator to prepare some large volume renderings using the aforementioned software volume renderer in yt.  In the past we've successfully made very, very large image renderings using yt -- Sam's even made one at 8192^2, although at extremely high resolution like that sometimes the lack of fidelity in the underlying volume renderer shows up; sometimes even artifacts in the AMR grid boundaries, but that's less common.  Making the very large volume renderings isn't too bad -- it scales roughly with the number of pixels, but we can dispatch many frames to be rendered at once on a cluster.  

There are a couple other, more important things to consider when making the big volume renderings.  For starters, the entire structure of volume rendering in yt was not really created to generate a series of images -- only a single image.  The idea was that you would prepare a specific image, make it, and move on.  However, for this project, I want to do a zoomin, or possibly a more complicated camera path.

Additionally, one of the first things that we did with the volume rendering was silly: we applied no normalization to the output images.  That was a mistake, I see now.  Part of the reason for this was uncertainty in the correct normalization -- the bias that the user wanted to apply may not be the natural bias from the image.  But more than that, because the rendering algorithm itself was some what holistically settled upon (the original implementation, which we used for shell-style renderings, was not a "correct" implementation of alpha blending) a natural mechanism for scaling did not immediately present itself.  One likely exists, possibly dependent on the field of view, I simply do not yet know it.  This will have to be rectified, because the mechanism used for scaling a set of images will have to be different than the mechanism for scaling an image in isolation, or else frames will jump in brightness during the movie's course.

The final thing that I wanted to change was to add support for stereo rendering.  Rather than repeat any of the amazing discussion from Paul Bourke's website, I'll simply direct you there.  Everything you ever wanted to know about stereo rendering.  (When I was a first year grad student, we actually bought a copy of his site to use locally -- it was our way of showing support for him putting it online, and it also came with a bunch of source code for example applications.)  I first attempted to apply the correct method for stereo, where the view direction is parallel and the total view frustum is shifted.

This did not work.  In fact, it made me realize that all this time, the yt method for volume rendering is in fact ... not really a volume rendering method inasmuch as it is a planar ray-casting method.  Typically when doing volume rendering, there's a perspective applied to the image: the rays all emanate from a single place, creating a frustum.  But for yt, we actually set up a single plane of vectors at the back of the volume and advance that forward across the image.  This is good and bad; it's good in the sense that it's more clear precisely what is going on.  But it's bad in the sense that correct stereo is more difficult.  (Of course, on Bourke's page he has a workaround that may work for this, but I have not yet attempted it.)  Here's a rough depiction of the different between the two methods.

Renderingmechanisminyt

The upshot is that stereo doesn't seem to work unless you go with the "toe-in" method that can cause eyestrain after a long time and shows visible parallax at the edges.  I'm not sure if this is going to be a problem, but because I am not right now eager to rewrite the rendering backend, this is the way it is for the moment.

To set up the stereo rendering, I separated out the rendering mechanism from the objects to be rendered.  Previously, there was a single VolumeRendering object that you could create, raycast through, then discard.  I created a new camera object that accepted a homogenized volume and would call "traverse" on that volume, feeding a back and a front point.  The Volume is then responsible for passing off fixed-resolution grids to the camera, which accumulates an image buffer by calling the ray traversal functions.  The front and back points are essentially the only thing needed to know this order, but the camera also stores its three orientation vectors and its position that describe it in 3D space.  By separating out these two conceptual objects, we undo some of the "single, carefully constructed image" bias that was in the original volume renderer.  (And, we open ourselves up to being able to use the Camera with a hardware volume renderer, should that day ever come.)

So now we have a camera, and it makes images like this:

C_0001

It's a little dim, but that's a task for another day.  The next step is taking that perspective and turning it into a set of stereo images.  To do that, I added a new class called StereoPairCamera.  It accepts a Camera object and turns it into two camera objects, where the final interocular distance is calculated relative to the image plane width.  As I mentioned above, this only operates via toe-in stereo, so it does this in the simplest manner possible: it moves each of the Left and Right cameras by half of the interocular distance away from the original location, and then recalculates a normal vector to point back at the original center.  Now we can generate left and right images:
Unfortunately, on my laptop (which is my primary/exclusive work computer) I don't have the ability to view these pairs.  To get around that, I wrote a simple stereo pair image viewer in OpenGL and imposed upon my friends that do have a stereo viz wall to test it out -- and after some fiddling with the interocular distance, we got what appeared to be workable stereo pairs.

The full code for generating camera paths as well as stereo pairs is already in yt, but the documentation is still being written; I might also clean up the interface a bit.  Additionally, at some point in the future, the issue of toe-in stereo versus correct parallel-frustum stereo will need to be dealt with; the last thing I really want to do is force people to only use a bad method for generating stereo pairs.  Hopefully that is something that can be dealt with at a later time.

Thanks to the wealth of resources out there for making this a relatively easy task: the aforementioned Paul Bourke website on stereo pairs, the PyOpenGL and PIL teams for making the image pair viewer easy, and everyone else whose work I've built on to make things like this.

L_0001

R_0001

Tagged python scipy viz yt

Some more fiddling with Chaco and AMR data

In a followup to a post from Friday, I fiddled a bit more with Chaco and the AMR data plotting today.  By mostly using components that I'd already created, I was able to make it switch between fields while panning and zooming.  By using Chaco and Enable a bit more cleverly, I managed to get a 25x zoomed in plot inset in the bottom left of the zoomed-out plot, and (as should come as little surprise) the necessary code was very short.  Here's what it looks like now, displaying a 512^3, refine everywhere, 7 level dataset:

L7_chaco

Panning and scanning and zooming through a dataset like this really is the only way to fly.  I'm now using it on data for my next paper, and I think it's about three lines of code from writing scripts to plot images based on clicking when I get to a scene I'd like.  There are a whole bunch of second-order fixes and improvements to make, but they have all been bumped to and unspecified time in the future now that it's functional to the extent that I can examine my data in detail!  But, I'll definitely push all of the changes necessary to get this type of interface back up into the main repository.  Perhaps somebody else will take it and run with it.

2D Data Visualization of AMR with Matplotlib and Chaco

With yt we have used matplotlib almost exclusively for 2D plotting of adaptive mesh refinement data.  This week, I implemented a plotter inside Chaco a plotting system designed for interactive exploration of data.

But first!  An interlude about how we visualize 2D AMR data.

Our data is naturally a patch-based format.  What that means is that higher resolution cells are grouped into grid patches -- this ultimately means that the data itself is relatively convenient to look at, and also for some post-processing computation we end up being able to make a lot of guesses at which data to look at or analyze based on its location in grid patches.  Additionally, these grid patches overlap at different resolution levels: finer resolution grids are inside coarser grids.

When generating a 2D representation, we typically will take either an orthogonal slice or a "projection" through the data.  (The additional spatially-organized 2D representation we have in yt is the oblique slice, but it's something of a special case.)  Back when designing the projection and slice algorithms, we decided it was the most efficient to construct a set of "flattened" data: by projecting through the data, we end up with a set of data that represents the finest resolution data available along every sightline through the domain.  This is represented by five values: (x, y, dx, dy, z) where x and y are in the plane of the projected domain.  For slices, we get something similar, but this time it's only the finest data available at any orthogonal slice.

Projection

The alternatively mechanism would be to generate fixed-resolution arrays and fill them in; this gives and takes various things, but ultimately while we provide support for this type of function, and it's the only way to project inline with the simulation (slices can be done adaptively inline), we prefer the "adaptive" projections.  An additional advantage, one that appeals to me particularly (as my AMR runs go to 30 or so levels of refinement) is that we only have to project once with the adaptive projections and we can zoom and pan as much as we like.  More on zooming and panning later!

In matplotlib, there are two natural ways to present this data: by constructing a set of rectangular patches that this data represents and then filling them with the appropriate value, or by pre-pixelizing, using a bit of our knowledge about how AMR data works, and then displaying an image.  The advantage of the first is that we are able to provide only a single set of data to matplotlib, once, and then we can adjust the limits inside the context of matplotlib and describe the data exclusively in the context of the domain.  However, there's a lot of overhead involved -- often we end up with projections that, with just the five values alone, run into the hundreds of megs.  (For very large simulations, particularly.)  Generating a second layer of objects, one for every cell at the finest resolution in the image plane, is probably going to add too much overhead.  So instead we run the pixelizer!  (Of course, at some stage matplotlib also runs a pixelizer, to turn the multiple patches into a single output image.  But it does so with a little less knowledge than we have about the nature of the data, so the generality costs would cost us.)

The pixelizer is relatively simple: for every datapoint in the projection or slice, if it falls within the bounds of our image, deposit it into the appropriate cells with which it overlaps.  If it fully overlaps, deposit the entire thing; otherwise, only deposit a fractional amount.  This process is relatively fast, but not without cost.  The biggest cost is actually in discarding the points that do not fall within the region -- the iteration itself!  Because our data is periodic, we have to check multiple times for most datapoints whether or not they are within a given region.  This also costs us.  Overall, however, the process is straightforward and fast.  (Of course, I do realize that the entire projection process: from the initial grid patches, to summing over an axis, to the compositing and the pixelizing, is something OpenGL was designed to do!  But, we do it all in C.)

So now we have a mechanism that takes multi-resolution AMR data and turns it into an image, given only the data itself, a resolution, and a 2D bounding box.  Typically this is used by yt to generate matplotlib plots.  This week, though, I have been working a bit on a tile display wall routine -- dispatching rendering jobs to different nodes and then having them feed back some portion of the render job to a display tool that broadcasts it to multiple monitors.  (More on that process in a future blog post.)

However, one thing that was always really lacking was an intuitive mechanism for interacting with a plot: in years past, I'd hacked up a matplotlib and wx-based GUI that had some neat features -- a slider to set the width and re-centering abilities.  But it was clunky, and kind of slow, and I hated maintaining it.  The book "wxPython in Action" was quite good, and wx is definitely a very capable toolset, but I ended up using the "pubsub" mechanism in wx and I eventually got lost in the events and handlers and couldn't find my way out.

Enthought's Chaco toolkit, though, was designed for interactive use -- whipping up GUIs around structured data (using the "Traits" interface) is easy and straightforward.  A couple years ago, they added the ability to get 1D data via functions -- so the data could be recalculated on the fly, as the range needed changed.  This week and last, with the help of the head Chaco developer, I was able to extend this to work with image data, as well.

First, a brief description of the UI I was hoping to get: Chaco very easily presents a "google maps" style interface to image data, where you can click and drag to pan and scroll your mouse wheel to zoom in and out.  I wanted to have it set up such that I could do this with the AMR data -- as you zoomed in, it would re-pixelize such that it was displaying a smaller portion at higher resolution, much like with normal google maps.  Almost all of this was available in the standard Chaco distribution -- the pan and zoom tools, the displaying of image data, and the setup for inserting callbacks between different components.  Displaying a single fixed image in Chaco, with the GUI and everything, is on the order of a dozen lines, many of which are boilerplate.

By adding a FunctionImageData class (almost identical to the FunctionDataSource class, except extended to 2D and subclassing Image Data) I was able to mock up the re-generation of data: as the range of data displayed changed, the data handed back was generated.  By inserting a callback to update the range that the data mapped to (i.e., not only generating higher resolution data, but mapping it back to a smaller data region, thus ensuring it fully filled the image plot) the entire experience was complete.

The FunctionImageData requires a callback function that receives a "low" and "high" value that describe the new bounds of the visible data, but I also had to update the "index" -- to do all of this together, I wrote this helper class:

 class ImagePixelizerHelper(object):
     index = None
     def __init__(self, panner):
         self.panner = panner
  
     def __call__(self, low, high):
         b = self.panner.set_low_high(low, high)
         if self.index is not None:
             num_x_ticks = b.shape[0] + 1
             num_y_ticks = b.shape[1] + 1
             xs = mgrid[low[0]:high[0]:num_x_ticks*1j]
             ys = mgrid[low[1]:high[1]:num_y_ticks*1j]
             self.index.set_data( xs, ys )
         return b
 

(Note that even though I sung the praises of Traits earlier, I don't really use them here.  That will happen, it just hasn't yet.)  The "panner" here is something that updates a pixelized image when the bounds change.  What the check for index inside __call__ does is mock up the re-generation of the image indices, if an index is supplied.  This also transforms the initialization of the image plot into a two-strep process:

         pd = ArrayPlotData()
         plot = Plot(pd)
         self.pd = pd
         helper = ImagePixelizerHelper(self.panner)
         fid = FunctionImageData(func = helper)
         fid._data = self.panner.buffer
         self.fid = fid
         bounds = self.panner.bounds
         pd.set_data("imagedata", fid)
  
         img_plot = plot.img_plot("imagedata", colormap=jet,
                                  interpolation='nearest',
                                  xbounds=(0.0, 1.0),
                                  ybounds=(0.0, 1.0))[0]
         helper.index = img_plot.index
         self.helper = helper
  
         fid.data_range = plot.range2d
 

So here I set up the array plot data, I set up the plot based on it, I create the pixelizer helper, then I add a default value to the pixelizer.  This is a crucial step, because it's how the image plot knows what shape the data is.  (There are some other things about the image plot I'm not fully sure of just now, and this is one I now sort of understand -- the relationship between mappers, indexers and plots.)  I've then set the helper to have access to the image plot's index, and the data range is set to the plot's data range.  At this point the two know how to talk together.  Below this I add a few tools for zooming and panning, and a colorbar too

At this point, our plot works.  And it's simple and easy and, best of all, very interactive.  Every time the data range changes the pixelizer updates and passes back a new image: so we've got the best of both worlds.  We only have to store in memory a relatively small amount of data, and it goes to the highest resolution possible.  Here's an example, where I am zooming into a dataset of mine.

Amr_chaco2

(download)

It could go on for a while -- this dataset is really quite deep (this was three scrolls, but I could probably do fifteen or so) but I think this gets the point across.  I was able to zoom in just by scrolling with my mouse, and I can click and drag.  It's intuitive, it's very very fast, and I only had to write a tiny bit of code to get it to go.  This is what I want!

I think that maybe it's time to revisit having a GUI.  I'd forgotten how much fun it is to explore my data, and if I can really get one going in such a short amount of code, it might be very effective to do this rather than constantly having to tweak scripts and scp images.  And with the vast amounts of data I am now sort of drowning in, this will be crucial.  The dataset I originally decided to do this for is just immense; I was struggling at coming up with a way to see it all, to see all the tiny features and structures, because it's at a time in the universe when such things are relatively far apart, but also quite numerous on the length scales we're simulation.  Using this kind of an interface will let me discuss individual halos with the other members of the project, and we'll be able to really see and share and converse, rather than just looking at static images over email.

Another huge advantage that only just occurred to me is that the same code I used to write the backend for the tile display wall could be used for remote visualization in this same manner.  By repurposing that to feed directly back into Chaco, instead of grabbing the arrays myself, we might be able to avoid any copying of data (even the pre-projected values) onto a GUI-able host.  That could be very awesome, too.  The entire Chaco/Traits/ETS toolkit is very powerful and I think I see it being an extremely valuable part of our future plans for interactive analysis.

(edit: evidently I messed up the options for code highlighting and for including all the images in a single slideshow, and from the post-edit area I can't figure out how to fix them..  I guess I have a bit more to learn about posterous!  Sorry about that.)

Tagged python scipy viz yt

Pyglet on Cocoa

As an update to an earlier post, it looks like Pyglet now works on Snow Leopard.  I was able to run several of the examples, with some minor glitches, even using my own hand-rolled Python 2.6 installation.
Tagged python viz

First Steps with OpenCL

I managed to get PyOpenCL to work as a projection engine for some AMR data.  This image is of a simple PyOpenCL driver I wrote using Enthought's Traits as a GUI.  The idea is that you can change some parameters, change the routine it calls, hit recompile, and it'll pass all the bricks back in.  I'd like to extend this concept, but first I've got to ensure that my OpenCL kernels are doing the right thing!

Pyopencl_projection
The kernel shown above (and visible as the subroutine "proj" in my OpenCL repo) does a simple projection.  The image might even be a MIP.  I've since been able to make it ray cast, but with a bunch of artifacts I haven't been able to track down, including one that has stymied plotting all three channels simultaneously.  Here's an image of the red channel from the ray caster, plotted against a different colormap:

Hi_0_gist_heat

So it's sort of working, but the artifacts will need some work before they're all sorted out.  Ultimately, I'd really like to support OpenCL as a rendering engine, as well as the current software engine we have now.  I think having the TVTK frontend that dumps to a software or OpenCL backend for ray-casting the full volume would be ideal.  I also think there's really something to the idea of having a generic engine (the simplest form of which I included a screencap above) that just passes data through a shader and plots the results.  I'm not the person to write one, I don't think, but I've heard rumblings one might be coming out of a pretty fancy coder up in the bay.

Adaptive Mesh Refinement Volume Rendering with Python and Cython

Over the last six months or so, I've been working on implementing a software volume renderer for adaptive mesh refinement data -- specifically, for 3D visualization of astrophysical data.  This has been implemented in the yt toolkit and it's been developed in the open.  I'll give a bit of background, and how various other people have had their hands in this.  There are two things that make me really happy about this new development: the first is that I have an actual use case for it, and I've been using it quite a bit for that this week (rendering my data!) and the other is that a number of people in the development community have had their hands in it, from start to finish.

The volume renderer we ended up with isn't the best -- it still has artifacts, and I can definitely see where it could use improvement! -- and it's not the fastest, but it's something that several of us worked on and it's something we can use to make viz of our data.  I think it's probably better thought of as a way of preparing a visualization, not a way of engaging in visualization.

Background

I conduct cosmological star formation simulations.  To study the formation of stars on cosmological scales, we use a technique called adaptive mesh refinement, wherein regions of interest are calculated at higher resolution than comparatively boring regions.  While this enables us to do awesome things like zoom in by a factor of over 10^12 in length scales, it also means handling the data is more difficult.  A couple years ago, my thesis advisor suggested I start work on an analysis toolkit written in Python.  So, I did.  It handles objects as described by physical regions (as opposed to code regions) and it abstracts away things like addressing of objects, data IO, generating fields, and so on and so forth.  This way, we could address things as "slice here through the domain" or "get me a sphere" and the toolkit would go on its merry way and take all the necessary steps to get the data and then get out of the way.

Over time, new features were built on top -- the viz capabilities mostly stayed where they had been since 2006, with a few additions but almost all were reductions to two-dimensions.  The only 3D visualizations we had were based on VTK, utilizing their HierarchicalBoxDataSet object.  The image below shows the yt "next-generation" GUI applying marching cubes to an AMR dataset, using the VTK operators and data structures.  There's a corresponding user interface, but it's still rough around the edges so I haven't included it here.

Star_vtk2

Late last summer, a colleague (and friend) of mine, John, implemented a mechanism that took the oblique-slice data object I wrote and made it into a "fixed resolution" oblique-slice.  By stacking up a bunch of these fixed resolution oblique slices, the equations of radiative transfer could be solved and a "volume rendering" could be obtained from it.

Casting A Couple Rays

John suggested that it could go a bit further, and so I started looking at actually performing direct ray casting through the simulation.  This is where a ray is thrown, and then it accumulates a value -- typically this is via a transfer function, where at every point the ray picks up emission and then via absorption it is diminished.  I'd written a single-ray caster back in 2008, so extending that to serve the new multiple-ray purpose wasn't terribly difficult.  However, what turned out to be the biggest difficulty was actually in the process of traversing grids.  In the adaptive mesh refinement implementation we utilize, some regions are covered by multiple grids -- "child masks" are applied to coarse grids to indicate where finer (better) data is available.  For the longest time, I had this mental block in my head that what I'd have to do is figure out a way to have a ray enter and exit a grid while traversing that very grid; I couldn't wrap my head around how to avoid this problem, or even how to handle it if I decided not to.  All I could see was that we had a bunch of rays.  These rays had to pass through every cell they intersected with, and they had to do it in order.

Raycasting

At my old institution, the wickedly smart Ralf K was the visualizer on-hand.  Besides being a code demon, he also wrote his thesis on the visualization of adaptive mesh refinement data.  During a conversation with him about this, he explained that typically what's done is that the grids are broken up into "bricks" -- non-overlapping, finest-data objects that fully tile the domain to be rendered.  (Ralf had also taken me under his wing and shown me techniques for visualization last Spring, before I graduated.  We prepared some visualizations of a simulation, and without exception his blew mine out of the water, even using the same software!)

The first attempt at a grid homogenization algorithm that I came up with was terrible.  It generated 64 bricks for every grid (It is not coincidental that this is (2^3)^2, I think) and it was slow.  Plus, with that many small grids, the advantages of viewing the cells as collections instead of individuals were basically nil.  This is where the project languished for a couple months.  We could make okay renderings, not great.  We would do the stupidest thing and linearly interpolate at N points along the ray between intersection points of cells -- even if it was only a glancing blow, the ray would interpolate N times between enter & exit.  This particular behavior hasn't changed, but it might some day.  There were artifacts, they were slow, and nobody really cared about making volume renderings -- for those reasons!

Homogenization

Over Christmas break, I was waiting for simulations to run, it was snowing, and I needed a fun project to spend some time on.  There are some papers that discuss it, but I am embarrassed to admit that most of them took it as written that the process of homogenizing adaptive grids is a solved problem.  I forget when it was, but I remember finally got over my inhibitions: I have had an aversion for the last couple years to writing code that I know is ugly, or that I know is inefficient, when I feel that things could be done better.  But, after I'd had enough, I just decided to let it go and I pushed forth.

I'm still not sure it's the right method.  But what I ended up doing was writing classes for GridFaces and ProtoPrisms.  The protoprisms sweep over directions until they encounter a grid face with which they intersect, and then they sweep in the next direction until the same.  Once a protoprism without any intersecting grid faces is identified, that prism becomes a "brick."  (The code, in yt/_amr_utils/VolumeIntegrator.pyx, is probably clearer than I am being right now.)

The final implementation I came up with only generates 2-8x as many bricks as grids.  It's sub-optimal, as it doesn't perform any minimum bounding box calculations, but it met my standards for a "success"!

Transfer Functions, Color Maps and Artifacts

So now we had a ray caster that worked.  All it needed was data, a transfer function to describe how to update the red, green, blue and alpha channels, and the will to push "go."  Also, the resultant images had artifacts, which was bad, but staring at the source code I couldn't see the sources.  That aside, however, making a transfer function was simply not easy, and I wasn't good at it.  So, using Enthought's TraitsUI, I sketched up a quick transfer function widget.  You can add new Gaussians, adjust their properties, and it calls a callback routine whenever "Run Routine" is clicked.  I was pretty happy with it, and it let me make some interesting images.  Seen here is a turbulent present-day star forming region, courtesy of my officemate David.

Tf_widget

The problem with this, however, is that I'm not really very good at generating color maps or schemes!  Fortunately, Jeff wanted to make some visualizations of data he was working with.  He came up with the idea of taking existing colormaps, in particular those provided by matplotlib, and using those as input to the transfer function.  The widths and the alphas were free parameters, but using the RGB combinations was a great idea.  Sam ran with this a bit, and I added some functions to add layers, and it all came together.

Sam, though, gets a substantial amount of credit for taking a look at the artifact problem.  He realized that rays were 'bunching up' at the edges of grids, so any isocontour (really, any non-zero emissivity, but we mainly use gaussians for isocontours) would show up grid lines that it crossed, very strongly.  The mistake was not only mine, but it was a mistake I came very close to identifying a few weeks before Sam fixed it.  (Note that I said: came very close to identifying, not fixing!)  With this fix in place, and Jeff's colormap ideas, we had a working volume renderer that made images we could all be pretty proud of.

Enstrophy

This image, made by Devin, is one that I think really shows off what makes me so proud of this volume renderer.  First, he defined a field that required extra cells to create a stencil; yt then fills the necessary cells from neighboring grids.  His function then creates a field.  (In this case, enstrophy.)  This field is then interpolated to vertex-centered data, homogenized and passed to the ray caster, where he then dumps back to disk and makes a movie.  I think it kind of looks like a Sentinel, from the Matrix.

Optimizations

The ray caster itself, along with the grid homogenization is written in Cython.  yt is a python package, but it includes hand-tuned HDF5 readers in C, and a number of the compute-intense routines are either in C, Cython or Fortran.  The usage of Cython for the ray casting has been an amazing boon.  "cython -a" is outstanding, as it drops out an html file with highlighted hotspots.

The zeroth pass of optimization was to try to make the algorithm as good as possible!  The first pass of optimization was to get rid of all of the yellow lines.  This wasn't terribly hard -- we mostly use pointers and C data types.  The final pass, though, was to utilize the Shark utility on OSX.  This is a free tool, provided by Apple, that can profile C/Fortran/whatever code.  Cython helpfully comments all the generated C code with the original line numbers, and it's relatively easy to figure out the correspondence between the Cython code and hotspots in the C code.
It turns out that the second order optimization actually revealed things I missed in the zeroth order optimization -- I was doing stupid and unnecessary things to interpolate values inside a cell.  After some tweaking, touchups and so on, the volume renderer actually sped up for some people by between 2x and 9x!  Much of this was due to culling cells that weren't inside the transfer function boundaries, removing calls to expensive functions like "fmod" and so on.  By multiplying by pre-calculated numbers, rather than dividing on-the-fly, we squeezed out maybe 10% or so.

At this point, I was happy with where the volume renderer was, and I was ready to start making movies and images.  So I did.

Seed008_blossom

This is just a still frame from one of my simulations; I have deliberately placed the isocontours here.  In the gallery (linked below) I have a zoom-in movie from a different simulation where the contours have been placed at evenly-spaced intervals.

Images and Movies

Things have progressed a long way!  Sam, Devin, Britton and I have all made some movies we can be proud of.  Sam has set up a gallery of his own and we've got an official yt vimeo channel and a gallery from that, too.

Next Steps

The very next step is going to be allowing the alpha value depend on a different field than the RGB values depend on.  John has suggested we give this a go to try to generate simulated images; I know this is Jeff's biggest priority, too.  I believe that this addition should be relatively simply, with a great deal of scientific payoff.

For now, though, I'm really, really enjoying using this on my own data.  That's why I wanted it in the first place -- to understand my data better, to reveal things about my simulations, to inform my intuition and assist my analysis.

There's a great body of research on the subject of volumetric visualization of adaptive mesh refinement data, by people that have thought about this a lot harder and in more depth than our group has -- and their results definitely underscore that!  Be sure to check out the amazing renderings available, and remember that the results above are still, one-off images -- not interactive explorations.  Those qualifiers aside, I'm still really proud of what we can make with this volume renderer, and I'm even more proud that it's all been a group effort, with the code all open source and available.

This writeup hasn't ended up being quite as technical as I originally thought, but I hope it's been enjoyable for you.  Be sure to check out the links to the galleries -- there've been some really great movies and images, showcasing some really fun scientific research.  Please also feel free to leave comments or questions -- and thanks for making it this far.

Essays Inspired by Microsoft’s Jim Gray, Who Saw Science Paradigm Shift

Dr. Wing has argued that ideas like recursion, parallelism and abstraction taken from computer science will redefine modern science. Implicit in the idea of a fourth paradigm is the ability, and the need, to share data. In sciences like physics and astronomy, the instruments are so expensive that data must be shared. Now the data explosion and the falling cost of computing and communications are creating pressure to share all scientific data.

There's a bit in here about the Worldwide Telescope and Google Sky, as well, along with discussion of Pan-STARRS and other new surveys.

The challenge for observational astrophysics may be more approchable than that for simulation-based astrophysics. A single simulation can generate terabytes of data, but the mechanisms for reducing and providing that data are still being developed -- although the Virtual Observatory folks have been taking this challenge head-on. Perhaps it's my own myopia, but when I look at different sets of simulation data, I see not only the astrophysical phenomena, but the different physics modules and the different methods for simulating them; when I look at the sky I see a coordinate system overlaid on top of a rich backdrop of astrophysical phenomena. The de facto organizational scheme presents itself.

Anyway, it's all a very interesting problem, but I think for now I will place my trust in the clever fellows with the organizational foresight, and I'll go back to making data that can be filed away nicely some day.

(Also, I think this might be the first time I've been interested in buying a book by Microsoft Press.)

Volume Rendering Widget

Media_httpfarm3staticflickrcom273541897493744dc637347fjpg_ubnfedyjxifhtdq

I wrote this widget this evening to edit transfer functions and toss a simulation through them. It's based on the TraitsUI system, which made it super easy to write, and the backend is a software ray caster I wrote. The data is courtesy Dave (my data is up next, but his is good to prototype with) and the rendering (512x512) is in the upper right.

The transfer function has several peaks, but I'm still tweaking how to represent alpha, so they don't all show up very clearly in the image in the upper left.

The UI is still a bit busy, but it's relatively straightforward, and I think it will help with editing and getting a good transfer function. The idea is to prototype it here and then dispatch the full rendering to a batch queue or some off-screen process.

Hollywood gives biologists a helping hand : Nature News

E. coli cytoplasm simulationThe movement of proteins through the crowded cytoplasm of a bacterium has been simulated.Adrian Elcock

Computer programs like those used in animated movies such as Shrek could soon be helping more cell biologists explain hypotheses — or even to make new discoveries, according to scientists presenting work in San Diego this month at the meeting for the American Society of Cell Biology.

"We want to be able to make predictions," says Adrian Elcock of the University of Iowa in Iowa City. "At the very least we want our models to reproduce known behaviours." Elcock is simulating the movement of proteins and other big molecules inside virtual bacterial cells. He built models from known data — including the atomic structures of proteins and concentrations of the 50 most abundant macromolecules in Escherichia coli — and then factored in how the molecular structure of each might cause proteins to stick to each other. His model nicely reproduces established data showing that green fluorescent protein diffuses approximately 10 times more slowly in the crowded environment of a bacterial cell than in a test tube.

Evidently Maya is now being used for biological visualization, and McGill even announced a "Molecular Maya Toolkit" to render proteins

Tagged science viz

NVIDIA Demos 3D Blu-Ray On 3D Vision - HotHardware

NVIDIA has been demonstrating a complete 3D movie solution to movie studios, press, and customers, consisting of a PC equipped with a GeForce® GPU and NVIDIA® 3D Vision™ active-shutter 3D glasses, as well as new 1080p, 3D LCD displays from Acer to showcase how consumers will experience this new 3D Blu-ray content once it is commercially available. In anticipation of an official announcement by the Blu-ray Disc Association and expected product unveilings at the upcoming Consumer Electronics Show (CES) this January, NVIDIA has been working closely with the world’s leading movie playback software developers, including Arcsoft, Corel, Cyberlink, and Sonic, to ensure seamless support for 3D Blu-ray titles when they are ready to ship in 2010.

Interesting -- does this mean a portable format for 3D scientific movies and a reasonable mechanism for displaying them in 3D on workstations is also coming? Right now the solutions are not very attractive, but perhaps this will help jumpstart that.

Tagged 3d viz