Outline of our program

This post is a simple step by step description of the main points in our program.

Setup scene
Build BVH over the objects
for #number of passes
  Trace the initial beams from the lightsources
  Calculate scattering events 
  Build BVH over the beams
  for #vertical pixels
    for #horizontal pixels
     Check for intersection with object BVH
      Collect intersected beams from the beam BVH
      calculate light from material
      calculate light contribution from beams

Scattering implementation

After implementing scattering in the participating media we had some problems when rendering highly scattering materials. Our first implementation used recursion, so with highly scattering materials we had a recursion depth of 100+. This caused a stack overflow in out program, which we solved by keeping a vector of beams we haven't checked for scattering yet, and then just looping trough that vector. The created scatterbeams are then added to the vector instead of recursively scattering.

void

Scene::addBeam(Beam beam) {

beamsToBeAdded.push_back(beam);

}

void

Scene::addBeams(){

int startSize = beamsToBeAdded.size();

//we run through the vector of beams that haven't been checked for scattering yet

while(beamsToBeAdded.size()>0){

Beam beam = (Beam)beamsToBeAdded[0];

beamsToBeAdded.erase(beamsToBeAdded.begin());

//we only scatter if the participating media is inhomogeneous or the      scattering value is greater than zero

//There is no reason to add the beams to the tree if the scattering coefficient is zero

if(hetero || beam.pm->getScatter(Vector3(0), beam.color) != 0){

if(beam.scatterDepth < maxScatterDepth) scatterBeam(beam);

//either subdivide the beam or add the current beam to the bvh.

if(useSubdivide){

subdivideBeam(beam);

}else{

beamvector.push_back(beam);

}

}

}

}

Another problem we encountered related to scattering was when combining scattering with the kd-tree for inhomogeneous participating media. When splitting the beams according to the kd-tree we had to make sure that only onescattering event occured across all the created split beams. This required some additional book keeping where each splitbeam knows the previous one, and can ask that beam if it's allowed to scatter.

Light contribution from each photon beam

The light contribution from a photon beam can be described by the following formula:

In our implementation we have implemented the basic idea in the following code segment:


for(int i = 0; i < beamHits.size(); i++)
{
BeamHitInfo bhi = ((BeamHitInfo)(beamHits)[i]);
//make sure the beams doesn't cross the object boundaries
if(ray.pm == bhi.beam.pm){
//not quite sure if this is the correct way to make the blur kernel?
float kr = (1/(2*PI * bhi.beam.radius));

Vector3 flux = (bhi.beam.flux/amountOfBeams);

//only use the color of the beam
if(useColors){
//multiply with 3 because each beam gets a third of its flux removed
flux *= 3;
if(bhi.beam.color == 0) flux *= Vector3(1,0,0);
if(bhi.beam.color == 1) flux *= Vector3(0,1,0);
if(bhi.beam.color == 2) flux *= Vector3(0,0,1);
}
int rayPoints = ray.getPointsUntilT(bhi.rayT);
float extMax = bhi.beam.pm->getExtinctionMax(bhi.beam.color);

float extMaxBeam = bhi.beam.pm->getExtinctionMax(bhi.beam.color);
float scatterMaxBeam = bhi.beam.pm->getScatter(bhi.P, bhi.beam.color);
if(useKDtree){
extMaxBeam = bhi.beam.maxExt;
}

//the flux reduction from the ray passing through the media
float trRay = exp((float)-rayPoints * extMax); 
//the flux reduction from the beam passing through the media
float trBeam = exp((float)-(bhi.beam.getPointsUntilT(bhi.beamT)) * extMaxBeam); 

//using the phase function to determine the probability that the light is scattered towards the ray direction
float angle = acos(dot(-ray.d, bhi.beam.direction));
float phase = bhi.beam.pm->getPhaseProbability(bhi.beam.color, angle);
float sinRayBeam = sin(acos(dot(-ray.d, bhi.beam.direction)));

//from the progressive photon beams paper
//lm = kr(u) * scatter(xw) * flux * Tr(ray) * Tr(beam) * (phase function(ray dot beam) / sin(ray, beam))
Vector3 color =  kr * scatterMaxBeam * flux * trRay * trBeam * (phase / sinRayBeam) ;
shade += color;
}
}

More heterogene media images

So basically we created a function, which would build a "wall" of 7x7 cells with different random scatter coefficients. The setup looks like this:






And the image after 68.000 photon beams:

And a 10x10 grid with ~550.000 beams:

Partitioning the media

We have been working on partitioning the media in a KD-tree to speed up the image rendering time as well as the contruction time of the photon beam map.

We define the space of the media with a bounding volume. This volume as well as any nodes/leaves in the tree are axis-alligned.
At the moment we find out wich axis to split along by just taking the longest side of the bounding volume. This should probably be changed, so we analyze the best split along all three axis instead of just one.
When the axis is found, just like the paper, we construct a graph describing the "thickness" of the media along the given split axis.


 

Say the split axis is the x-axis, see image above. To construct this graph, we move along the x-side of the bounding volume of the tree-node in small steps.
Each step is an x-coordinate and this defines a plane or "slice" through the media perpendicular to the x-axis. We take a number of samples in this slice, the red dots in the image above. We keep track of the largest (max) and smallest (min) extinction value found in each slice, and the value
graph(x-coordinate) is then max-min.
The graph could the e.g. look something like this:



Meaning the media is thicker with higher x-values.

To find the x-coordinate to split along, we attempt to solve the largest empty square problem.
The idea is that we want to create a split, which would divide the media in two: a volume with high extinction and a volume with low. We want to find the largest square above the graph, as shown in the picture.
A division of the above graph could e.g. be ~1/3 in, as marked with the blue line. This will give us the x-coordinate and plane to use to subdivide a node in the kd-tree.
To determine whether to subdivide or not, we have a cost function:
   area of rectangle - 1
If this is > 0, we split.

We now need to use this tree when rendering and constructing the KD-tree, to see if it speeds things up.

Speeding up inhomogeneous participating media

The second part of out study group is implementing a paper to speed up the rendering of inhomogeneous participating media. The paper can be found here:

http://nis-lab.is.s.u-tokyo.ac.jp/~egaku/sigasia10/abstsigasia10.html
or
http://dl.acm.org/citation.cfm?id=1866199

The basic method is to partition the participating media into a kd-tree based on the extinction coefficient, and one of the key aspects of the paper is a scheme to find a suitable partition.
Our plan at the moment is to first add the functionality to render the participating media with a kd-tree, and then implement the partition scheme. We will probably just start with a uniform partition to get the rendering working correctly, and then afterwards try to get speed improvements with the good partitioning.
The reason inhomogeneous materials can be slow to render is because of the woodcock tracking used to get  the amount of color absorbed in the material.
Woodcock tracking works by taking smalls steps bases on the maximum extinction coefficient in the material, and then use the extinction coefficient at the current point to determine if a scattering or absorption event occurs. The problem with the method is that if the extinction coefficient varies a lot over the material, a lot of unnecesarry steps can be taken. So the idea is to partition the medium into partitions with similar coefficients and then takes steps according to the maximum extinction coefficient of each particular partition. 

woodcock tracking

Today I have tried to compare our rendering of homogeneous media with the same material settings using woodcock tracking and the other calculating the extinction directly. They do (luckily enough) look quite similar, but their was a huge speed difference between the renders. First a render with our homogeneous system of merlot in low concentration:

 This is a render with just one pass, but because we can calculate the amount of light absorbed in the material with e^(-t*extinction) we can get a very clear result with just one pass. The colored dots that can be seen is the result of rays being reflected on the surface of the liquid, which we are doing with russian roulette which requires more samples per pixel to get a smooth image. The woodcock render of merlot in low concentration looks like this:

This is after 35 passes of the render. This render is done using woodcock tracking, and it only approximates the correct light absorbtion in the material. It is an unbiased method, so given enough passes it should convert to a correct solution, but because of the randomness of the algorithm is requires a lot more passes to convert to a clear result. It does however look quite similar to the homogeneous render, so our next project is to speed up the rendering.

homogeneous coke 1 pass:

woodcock tracking, coke 68 passes: