"clipka" <nomail@nomail> wrote:
> http://graphics.stanford.edu/papers/bssrdf/

I just worked through a good share of this the other day.  You'd sure be in luck
if I had a better grasp on it...


> Note that the paper doesn't specify a formula for "alpha_prime", which is
> supposed to be the "reduced albedo". In most sources I found this defined as
> "sigma_prime_s / sigma_prime_t". Some source, however, specifies it as
> "sigma_prime_s / sigma_a", which makes the code yield *similar* results as the
> approximation does with the "standard" definition of alpha_prime.

You can find an implementation (PBRT) here.  According to their code,

Spectrum ap = (sigma_sp / sigma_tp);

And they get good results...


> So, here's my C++ implementation - anyone able to spot the difference?

I'll try my best.


>     double sigma_prime_t = sigma_prime_s + sigma_a;
>     double F_dr = (-1.440 / Sqr(mu) + 0.710 / mu + 0.668 + 0.0636 * mu);

eta?  Fair enough.  Where did you get the sigma_prime_s and sigma_a from?  Does
this require evaluation of the mean cosine, or do you assume isotropy?  Should
these be vector parameters for the rgb values?

>     double Aconst = ((1 + F_dr) / (1 - F_dr));

These are all promoted correctly to double precision, but I try to avoid getting
lazy here since I've done some pretty stupid things in the past...

>     double dist;
>     double Rd;
>     double phi_in, phi_out;
>
>     VDist(dist, in.IPoint, out.IPoint);
>     dist *= sceneData->mmPerUnit;
>
> #if 1
>     // full BSSRDF model
>
>     // TODO FIXME - why doesn't this give roughly similar results as the BRDF
> approximation below??
>
>     double Dconst = 1 / (3 * sigma_prime_t);
>     double alpha_prime = sigma_prime_s / sigma_a;  // [1]

See note above regarding choices.

>     double sigma_t_r = sqrt(3 * sigma_a * sigma_prime_t); // [2]
>
>     double  zr = 1.0 / sigma_prime_t;
>     double  dr = sqrt(Sqr(zr) + Sqr(dist));

"dr = ||x-xr|| is the distance from x to the real source, and dv = ||x-xv|| is
the distance from x to the virtual source."  I need to take a closer look at
this section.  You're sure of this formula?

>     double  zv = zr + 4 * Aconst * Dconst;
>     double  dv = sqrt (Sqr(zv) + Sqr(dist));

Same as above.

>     // calculate Rd
>     double m  = alpha_prime / (4 * M_PI);
>     double t1 = zr * (sigma_t_r + 1.0/dr);  // [3]
>     t1        *= exp(-sigma_t_r * dr);
>     t1        /= (sigma_prime_t * Sqr(dr));

For a second I thought there was an error in your simplification, but except for
the zr, this looks good.  This does seem much more symmetric, so whether or not
it works, this form would seem more logical.

>     double t2 = zv * (sigma_t_r + 1.0/dv);
>     t2        *= exp(-sigma_t_r * dv);
>     t2        /= (sigma_prime_t * Sqr(dv));
>
>     Rd = m * (t1 + t2);

The code that you've provided looks good for the most part, but I'm assuming
that you have correct input and output.  This doesn't look easy to debug since
there aren't many sanity checks for intermediate values.

 - Ricky

Post a reply to this message