Modified this post 2012/6/21
I have posted images of particles. It might be interesting to use this one.
"Noise Reduction for Path Traced Imaging of Participating Media" (http://www.eurasip.org/Proceedings/Eusipco/Eusipco2011/papers/1569426055.pdf).
Another thing I would like to mention is about "Accelerated Light Propagation Through Participating Media"(http://isg.cs.tcd.ie/cosulliv/Pubs/LeeLightProp.pdf).
I found a little mistake in updating pos. I think it has to be
pos = pos + dir x ( -log(u) / (sigma_a+sigma_s) ).
The authors terminate a ray if ta < ts, which I thought was interesting. The parameters ta and ts are given as
ta = -log(u1) / sigma_a
ts = -log(u2) / sigma_s.
However, a little calculation shows that this is exactly the same as performing Russian roulette with sigma_s / (sigma_a + sigma_s) so there is no need to compute log() every time.