Search This Blog

Sunday, March 27, 2011


Anybody a CMake guru that can have me rapid-fire dumb questions their way? I've got ideas on how I'm going to make my API work, but now I need to really nail down the build system. I have a number of external dependencies that I can either require the user to download and install on their own, or I can find some way of grabbing all of them for the user, and putting it together into a package. The ideal method would involve my telling CMake what those external dependencies are and where to find them, and let it handle the chore of making sure that they are properly installed on whatever system the user has, offer to download and install them if they aren't available, or fail gracefully if they aren't. Note that I really don't want to be responsible for updates to external dependencies; as far as is possible, I want them to 'just work'. That will let me concentrate on coding, rather than packaging.

Tuesday, March 22, 2011

Current plan

After reading through various papers, and talking to Dr. Jaggard, I've realized that I can't tackle everything there is in the world of optics at one time.  I have to pick and choose my fights.  For that reason, I'm going to try to implement the following capabilities in the following order:
  1. Unpolarized light, in a world of perfect insulators and perfect conductors.  This is probably the simplest physical model of light, and is due to the fact that light is just  Electromagnetic radiation that we can see.  Since it has an electric field, that means that as soon as a wave hits a superconducting surface, it is absorbed.  That means that none of the objects in this model will have any transparency what-so-ever, except for the air, which I'm going to model as a perfect vacuum.
  2. Polarized light.  The next simplest thing to work on is polarized light; in this case, the most general model is elliptically polarized light, which can be analyzed using Jones Calculus.  The handy part about this is that if you squint a little, Jones Calculus starts to look like a rotation matrix that a first-year student has messed up somehow (it doesn't have to be skew symmetric).  Now, the reason why I need to get polarized light working quickly is because my primary interest in Fourier optics is for modeling radio propagation, where my antennae are all dipoles.  Two dipoles that both point in the same direction will be able to transmit to one another, but two dipoles that are 'crossed' will not be able to transmit much energy to each other (this is a simplified model, if you are a radio engineer reading this, don't shoot me! :)  OTOH, if you want to improve my model, feel free to comment here, or better yet, enter in an issue at the issues page.)  Note that as long as the world is full of perfect conductors or insulators, the materials in the world will not be able to polarize light; the light is created in a polarized state, and all we need to worry about is if the receiving antennae are in the right orientation to receive the incoming information.
  3. Hybrid raytracing/Fourier Optics model.  The problem with performing the Fourier transform across all of space is that it is expensive to do, with no guarantee that you'll really use all of the calculation that you so expensively performed.  In fact, once you're past about 10 wavelengths from where your photon packet is, you no longer really need to think about the wave nature of light; the ray cleanly missed the object.  One solution to that would be 'fat' rays; project a ray just as you would for ray tracing, but make each ray a solid rectilinear beam.  Do the FFT on the volume that the beam travels through, instead of all of space.  This will need to be optimized so that when the beams overlap in space, you don't waste time recomputing the FFT over that region, but since I know what the viewing frustrum is like, and how big my beams are, I can easily calculate the region of space near the camera where they will overlap (yes, you can develop edge cases where they don't overlap, but my suspicion is that the time wasted in determining if you overlap or not will be greater than just doing the FFT and throwing away anything close-by that you don't use)
  4. Imperfect conductors At long last, we arrive at a more realistic materials model!  The index of refraction can take on any value, because the permittivity and permeability of objects can vary.  Note that at this point, things get interesting (meaning hard), because both the permittivity and the permeability are rank 2 tensors.  This means that the permittivity and permeability of some small volume depends on where it is (a volume within a chunk of glass is different from a volume in air), and on the orientation of the volume (this affects how the material polarizes light that passes through it).  In addition, the permittivity and permeability affect how quickly light is extinguished within the material; if you think about what THAT means, it means that both must be frequency dependent, otherwise we'd never have colors.  Finally, you also have reflection and refraction.  Thus, modeling this is the most complex problem to date. 
The hard part with all of this is determining how to correctly model light, and to correctly model the permittivity and permeability of the material.  In addition, my need to perform ray-tracing directly affects where my code ends up in the OpenGL pipeline.  I can do all of point 1 in CUDA, and then emit the answer straight into a fragment buffer, and let it do all the heavy lifting of displaying the output.  However, the fragments won't have any information on how the rays went through space; I need have intersection information, as well as be able to generate new rays when there is reflection or refraction.  That suggests that I need to tap into the pipeline after the vertices are known, which means after the vertex shader stage (or geometry shader stage, if there is one) is done.  Taking everything together, my current plan is as follows:
  1. Define a world coordinate system from which all orientations are derived.  For the simple reason that I'm most familiar with it, I'm going to define a right-handed coordinate system where 'z' is 'up', 'x' is to the 'right', and 'y' is 'forwards'.  Along with this, I define the angle θ and ϕ, where θ ranges over [0,π] and 0 is parallel to the z-axis, and pointing 'up', and where φ ranges over [0,2π), and 0 is along the x-axis, pointing right.  ϕ = π/2 points along the y-axis, pointing 'forwards'.  
  2. Define every ray of light has 4 vector (intensity, frequency, θ, ϕ).  The downside of this model is that if a ray has multiple colors, then I will have to fire the same ray through the same medium multiple times.  Unfortunately, I can't think of a better method of handling this.  I will initially treat all rays as having an intensity of 1, a single frequency, and linearly polarized in the vertical direction.
  3. Model each cubic volume as having two rank 2 tensors, one for the permittivity, and one for the permeability.  Initially, these will be set to either the permittivity and permeability of free space, or to permittivity and permeability of perfect conductor.
I still have a lot more to think about, especially the ray-tracing part.  Too much work to do!

Talked to Dr. Jaggard

I talked to Dr. Jaggard yesterday about Fourier optics.  He gave me some suggestions, including looking at Classical Electrodynamics by J.D. Jackson, and the journals Radio Science and IEEE Transactions on Antennas and Propagation.  So this morning I checked out Classical Electrodynamics, and quickly read the chapter on propagation.  Unfortunately, I don't think that that chapter will be of much use; it breaks up the regions of radio propagation into very near, intermediate, and far, and then makes approximations for each of the regions so that you can solve them by hand.  I'm not solving anything by hand, so I don't need those approximations.  So at this point, I need to re-derive the equations.  This isn't particularly difficult as I'm going to be using many different books as guides, but it does increase my workload.

I'm also trying to figure out what translation and rotation in Fourier space really means.  Ideally, I'd find a general solution for an object, and then do an affine transform on its Fourier transform to get any translated, rotated, or scaled version of it that I want.  I don't know if this is mathematically possible, or it is better to do a pure FFT each time the environment changes.

I'm also trying to see if it is possible to perform limited Fourier transforms; that is, since I know the wavelengths I'm interested in (which is always some range), is it possible to reduce the amount of calculation necessary, so instead of integrating over all of space, I only integrate over those portions that are near the wavelength.  This may not be possible simply because of resonance; if I have a cavity (like the inside of a hallway), which happens to be a multiple of the wavelength, then there will be bright and dark spots.  If I don't include the higher spatial frequencies,  I risk making the hallways look 'blobby' rather than rectilinear, which means that resonance won't occur.

Most importantly, I still don't know if I can take one 3D spatial Fourier transform, and then convolve my cameras/lights whenever and wherever I want.  I think I can, but I can't state that for certain.

Monday, March 21, 2011

Design process

Designing a new API is one of the hardest things I have ever done.  Whatever you come up with, you're stuck with it forever... at least, once you hit 1.0.  Fortunately for me, I'm not anywhere near there yet, so I can make modifications as I see fit.  The downside is that I can see many, many different ways of designing the code, each of which has its own good and bad points.  Here are some of the problems I'm dealing with:
  • Choice of language - When choosing a language to write code in, you have to choose between performance and convenience, availability and modernity, etc.  Modern languages like OCaml and Python are powerful, and have sufficiently strong semantics that once you get your code compiling, you've probably solved most of your bugs already.  Unfortunately, their performance is pretty poor, and there is no guarantee that you'll find compilers for various platforms.  That forces your hand.  In my case, I had to choose between C and C++.  I've decided to use C, specifically C99, because it is widely supported, and because many of the higher level languages have the ability to bind to a C API.
  • Maintenance of state information - Next is how to maintain state information across function calls.  Older libraries maintained state information in an user-inaccessible global state variable or variables.  The problem with this approach is that more and more programs are multi-threaded; if each thread requires its own state information, you need to develop a method of doing so that is efficient, yet invisible to the user.  This becomes MUCH more problematic if a user wants different threads to share a context, or to have multiple contexts, or has multiple GPUs.  For that reason, my API has an opaque context pointer type, that points to a context that is allocated on the heap.  You can have an arbitrary number of contexts, and share them in any manner you wish.
  • Allocating & deallocating memory - malloc() and free() may be the C standards, but they aren't necessarily the most useful calls for an API, especially one that is designed to take advantage of a GPU.  Towards that end, I've designed my own memory allocation and deallocation routines.  They are designed to make debugging simpler, and to ensure that the necessary memory alignments to make use of the GPU are adhered to.  
  • Primitives - This is the part of the API that I'm having the most trouble with at the moment.  I have several different options, none of which are particularly good:
    • NURBS + material specification - This is by far the most general method I can think of.  Every object would be specified as a NURB surface, along with some information about the material it is made of (index of refraction, polarization information, transmittance and reflectance information for different optical frequencies, etc.).  With this information, the library performs the full simulation (Fourier transform, convolution, inverse transform) to determine what light does when it travels through the material.  There are a number of problems with this method though:
      • NURBS don't handle sharp corners well; you either have discontinuities, which can't be transformed, or you have very curved surfaces, but no sharp corners.
      • I'd have to figure out an analytic solution to a NURB surface.   While this can be done, I'd really rather not.
    • Library of primitive objects - I could define a small set of primitive objects, and require the end user to develop more complex objects from these primitives.  While this is conceptually simple (I get to choose the primitives, so I'll choose ones where the Fourier transform is easy), it means that the end user needs to build their objects from the library I provide.  The biggest problem with this is the interfaces between objects; regardless of how hard a user tries, there will be gaps between objects.  Those gaps will lead to numerical instabilities.  Worse, some objects require very large numbers of primitives; for example, how would you model a light bulb?  To a first order, it is a sphere within a sphere, but the neck is not a sphere, a cylinder, or anything else.  The only way to model it is via successive approximation, which can easily run out of memory.
    • Voxels + 3D FFT - A completely different approach from the analytical one is to use the discrete Fast Fourier Transform across a 3D voxel array.  The good part about this approach is that regardless of shape of the primitive, you can always find a discretized version of it to perform the FFT on.  Even better, if you tap into the OpenGL pipeline at the appropriate location, you can pull out a 3D texture, which means that the hardware did all the heavy lifting for you.  You also get all of OpenGL for free; you're tapping into the pipeline after it's been processed.  But on the downside, the smallest object you can handle is one voxel in size; the problem with this is that the wave nature of light is most important for objects that are on the order of the wavelength of light.  Visible wavelengths are on the order of 380nm to 780nm (, which, thanks to the Nyquist Theorem means that our voxels need to be on the order of 190nm on a side.  Assuming a laughably small 1 bit/voxel storage requirement, that means 1 cubic meter will require ~2^64 bytes of memory.  Clearly, a pure voxel approach won't work.
  • Rotation/translation - Beyond the problems outlined above, I don't know if the mathematics remain unchanged when an object is rotated or translated.  Thus, precalculating what a sphere looks like may be pointless simply because I can't move the sphere to a different location in space, and expect it to remain the same mathematically.
At the moment, I have no idea what will be the best method.  If anyone has any ideas, post in the comments!