Can someone check what's wrong with the following code? It's supposed to
implement the formula (4) from section 2,1 of this paper:

http://graphics.stanford.edu/papers/bssrdf/

The #if'ed-out section implements the approximation given by the first formula
of section 2.4 of the same paper, and seems to work perfectly.

I've been reviewing the code over and over, comparing it to the formula given in
the paper, another implementation, and various tutorials and papers further
exploring on the subject, but still can't make out why it doesn't seem to work.

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. However, I
can't imagine that the BSSRDF paper uses two concurring definitions, and anyway
the results are still a deal off the mark, and also seem to get closer or
further off depending on material parameters. (See [1] in my code)

Also note that formula (4) seems asymmetric regarding "zr" and "zv"; in fact,
various other sources add "zr" as an additional factor to the first term in the
angular brackets, making the formula symmetric. Strangely enough, this factor
does not seem to have any noticeable effect. (See [3] in my code)

One source also had a factor of 2 instead of 3 inside the root in "sigma_t_r".
Again, I did not see any noticeable impact. (See [2] in my code)


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


-----

    double sigma_prime_t = sigma_prime_s + sigma_a;
    double F_dr = (-1.440 / Sqr(mu) + 0.710 / mu + 0.668 + 0.0636 * mu);
    double Aconst = ((1 + F_dr) / (1 - F_dr));
    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]
    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));
    double  zv = zr + 4 * Aconst * Dconst;
    double  dv = sqrt (Sqr(zv) + Sqr(dist));

    // 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));
    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);
#else
    // uniform illumination BRDF approximation

    double  alpha_prime = sigma_prime_s / sigma_prime_t;

    // calculate Rd
    double  root = sqrt(3.0 * (1.0 - alpha_prime));
    Rd = (alpha_prime / 2.0) * (1 + exp((-4.0/3.0) * Aconst * root)) *
exp(-root);
#endif

    VDot(phi_in, *vIn, in.INormal);
    VDot(phi_out, *vOut, out.INormal);
    phi_in = acos(phi_in);
    phi_out = acos(phi_out);

    *sd  = 1 / M_PI;
    *sd *= ComputeFt(phi_in, mu);
    *sd *= Rd;
    *sd *= ComputeFt(phi_out, mu);

-----

Post a reply to this message