Hi everyone,

I've tried to update my ray tracer to handle Monte Carlo path tracing rather than standard ray tracing, and I have got it working but I'd like to check my algorithm over as I'm unsure about a few things. Here's what I'm doing:

Base rendering routine: For each pixel (loop through y, x) and for each sample (loop through sample number, e.g. 100 rays/pixel) shoot a ray through a random point in the pixel, add the results and then divide by the number of samples.

My trace function does the following:

- For each initial pixel ray, assign an attenuation vector with the value (1, 1, 1) which will be decreased as the ray propagates

- If there's no intersection, return a background colour

- If the intersected shape has an emissive colour, return that colour with 100% probability

- Define the following:

float kd = shape->material.diffuseColour.sum(), ks = shape->material.specularColour.sum(), kt = shape->material.transmitColour.sum();

...where '.sum()' adds together all the components of the vector

- Choose whether to spawn a diffuse, specular or transmitted ray like this:

diffuseAttenuation = attenuation*shape->material.diffuseColour; //Per-component multiplication

diffuseProbability = kd/(kd + ks + kt);

if (diffuseAttenuation.findmax() > cutoff && random(0, 1) < diffuseProbability) //'.findmax()' returns largest component of vector

{

- find a random (cosine weighted) direction

- spawn a recursive ray in this direction (its attenuation vector is now equal to diffuseAttenuation)

- multiply it by the current diffuse colour and return.

}

else

{

specularAttenuation = attenuation*shape->material.specularColour;

specularProbability = ks/(ks + kt);

if (specularAttenuation.findmax() > cutoff && random(0, 1) < specularProbability)

{

- spawn a recursive ray in a perfectly reflected direction (its attenuation vector is now equal to specularAttenuation)

- multiply it by the specular colour and return

}

else

{

transmitAttenuation = attenuation*shape->material.transmitColour;

if (transmitAttenuation.findmax() > cutoff)

{

- spawn a recursive ray in a perfectly refracted direction (its attenuation vector is now equal to transmitAttenuation)

- multiply it by the transmission colour and return

}

}

}

...OK, so here's what I want to ask about:

- I've seen some sites/papers say that I should divide each ray by its probability, e.g. the colour returned by a diffuse ray should be divided by diffuseProbability, etc... - should I be doing this?

- Are there any things you can see wrong with my algorithm? I don;t have any books on this so I'm just going based on what I can lift from papers and various PDFs from unviersities, none of which seem to tell the full story.

Thanks!

p.s. here's a picture it rendered which took about 45 minutes at 500 samples per pixel. Is this the kind of speed that I should be expecting with this type of path tracing?

http://i.imgur.com/hdVfmXh.png

:)

Seems more or less correct, it's the basic workflow of a monte carlo path tracer. Usually it's implemented iteratively though for performance. The diffuseProbability and specularProbability seem to be heuristics that are not physically based, but yes, it is the general idea - instead of launching a reflected and a refracted ray, tracing them, and multiplying their contribution by their respective probabilities (= exponential number of rays), you just select one at random according to their probabilities and give it full contribution, then they naturally average out to the correct weights.

That means that when you branch into a perfectly reflected ray, that code is already going to be running (and a new sample produced) with probability specularProbability, so the attenuation of the new ray should not be multiplied by the specularProbability (which you are correctly doing, I think). Mathematically you should also account for the distribution of the possible ray directions you are choosing from in that branch, but when reflecting specularly there is only one possible direction either way (unlike diffuse or glossy reflection, for instance).

- I've seen some sites/papers say that I should divide each ray by its probability, e.g. the colour returned by a diffuse ray should be divided by diffuseProbability, etc... - should I be doing this?

No, as I mentioned above. You are already implicitly dividing by the probability (in the limit) by branching on a random variable. But I believe they might be referring to "russian roulette" which is a technique used to reduce noise in the final image by giving more importance to the first few bounces and probabilistically killing rays depending on their intensity (attenuation), while keeping the result the same. For instance, if you have a ray that after a few bounces, has intensity 0.01, it isn't going to contribute as much as a ray that's still at 0.95 intensity. So what do you do? At each bounce, generate a random number, and if it's less than the ray's intensity, keep going, otherwise discard the ray. Mathematically it comes out the same in the limit, but the overall effect is that the ray tracer focuses on rays that are more likely to make a large contribution to the final image, which means less noise. In practice, though, you only enable this after a couple bounces, otherwise you end up focusing on the "too high contribution" rays and actually worsen the noise. But anyway. Most likely the paper was implementing something differently and needed to divide by the probability to account for that. It's impossible to tell.

By the way you might want to check those branches. In your code here, specular reflection/refraction occurs only if diffuse reflection is not selected. So that means that the probability of the ray reflecting specularly is actually (1 - diffuseProbability) * specularProbability, I believe. This may not be what you wanted.

- Are there any things you can see wrong with my algorithm? I don;t have any books on this so I'm just going based on what I can lift from papers and various PDFs from unviersities, none of which seem to tell the full story.

Welcome to the wonderful experience of finding source code samples (with explanations) for path tracers. Almost all the resources you will find are either wrong, misleading, incomplete, contradictory, unavailable, too complex, or even voluntarily obfuscated (yes, I'm looking at you, smallpt!). You really just have to sift through it and glean little bits of knowledge here and there and working out what makes sense and what doesn't. You will probably be in constant doubt throughout the whole process, wondering if you missed a constant here, a factor there, and realistically nobody will be able to help you efficiently at this level because your questions will be so intrinsically linked to your own understanding of path tracing (and to your own code) that they are more likely to confuse you than help but pretty pictures make it worth it! Trying to derive things "from the math" can also help a lot in clearing misunderstandings, but can also serve to confuse you even further, so tread carefully.

Well it looks good (at least it doesn't look obviously wrong, beyond that there's nothing else to say without a reference render). 45 minutes for 500 samples is a bit slow, but if your code is single-threaded and is not yet optimized it seems reasonable and not awfully slow. I would at least try to grab the obvious speedups (multithreading and the easy optimizations) so that you can prototype as fast as possible, it's definitely worth it.

Thanks for the detailed reply Bacterius :) So by iteratively I guess you mean rather than spawning a new ray with depth + 1 I would loop through all the depths and change the existing ray's direction/proeprties, all the while adding up the colour? That sounds like a good idea seeing as how I only spawn one ray from each intersection, although I'm using cumulative attenuation rather than a ray depth. Maybe I could do it in a while loop?

Yes, what I'm doing now is giving whichever ray is spawned full contribution to the sample.

Yeah, I'm a bit confused about my probabilities. Since I can only choose one secondary ray I have to check each case one after the other, so I can't see any other way of doing it (unless I limit objects to being only diffuse, specular or refractive).

It's currently running on a quad core 1.8Ghz laptop and I'm using OpenMP to multithread, but there aren't any other optimisations.

Thanks again :)

Alvaro - I looked at photon mapping too but it seems that it would need KD trees to implement and I'd rather not go there given that I'm still doubting myself about my basic implementation!

Ah right, I'll take a look I've seen enough papers and powerpoints that it will be nice to see an actual book!

In the meantime I've updated the path tracer to be simpler (each surface is now either diffuse, refractive, specular or a light) so that I can be sure it's rendering properly.

I also made it iterative, which actually makes it a lot simpler and intuitive (for me, anyway). Right now every ray starts off as (1, 1, 1) and it then enters a while loop based on the magnitude of the largest component of the ray's colour. Every time it hits an object the ray's colour is multiplied by the object's colour and it is given a new direction based on the BRDF of the surface. If it hits a light then the loop is broken and the ray colour is multiplied by the light colour, same for if it hits the sky. The picture it produces is basically the same (yes, I changed the colour of the sphere on purpose!!). This took 205 seconds at 50 samples per pixel.

I also now understand the Russian Roulette method for getting rid of certain rays, and I get why you divide by the termination probability. I'll try and see how it looks with textures too :)

500 samples now takes 30 mins to render:

[attachment=21415:weiut.png]

OK, I've been thinking some more about combining specular, diffuse and transmission. Rather than doing it using different colours for each and comparing their magnitudes, which is what I did above, I think this makes more sense:

- Each material has a single colour

- assign a probability to each of diffuse, specular or transmission for a given material which total to one, for example:

[material shiny_yellow]

colour = (0.8, 0.8, 0.2)

pd = 0.1

ps = 0.9

pt = 0

[material dull_glass]

colour = (0.9, 0.9, 0.9)

pd = 0.4

ps = 0.1

pt = 0.5

...and then in my code, do the following for an incoming ray:

- Find hit point, multiply ray colour by surface colour, get the three probability values

- generate p = random [0, 1] number

if (p < pd){ray is scattered diffusely}

else if (p < pd + ps){ray is scattered specularly}

else {ray is transmitted}

...in my head, this makes sense. I hope I'm not being stupid! Also, I still don't understand why some papers are telling me to divide by probability density functions left, right and centre.

That's a good point, I know that it's material specific but for simplicity I'm leaving that out for now! Also I think I finally understand the whole probability division thing - you need to do it if you are choosing whether or not to terminate a ray as opposed to which direction your current ray continues - is this right? So since I'm propagating all my rays until they run out of energy or hit a light, I'm not introducing a bias whereas if I was terminating them based on their radiance, I'd have to weight accordingly. Damn, it's so annoying trying to figure this stuff out from such contextless sources! Makes me wish there was a course on this I could take at university.