A few questions about the PBD implementation

[tomvbussel]

Hey,

I've been reading through the implementation of Photon Beam Diffusion (which is a really nice upgrade over the dipole in pbrt-v2 by the way), while comparing it to the paper and the Matlab reference implementation. I still have some questions about the code and I think I've found a few mistakes.

1. Is the normalization constant 1 / (4 * C_phi) missing in the BSSRDF? It seems that it was replaced by sigma_t which I can not find in the paper.
2. Are the 1 / pi term and the two F_t terms in the BSSRDF computed using the BSDF? You seem to both add a LambertianReflection (for the 1 / pi term, I guess) and a SpecularTransmission (for the F_t term), but when evaluating the BSDF this results in (1 / Pi + F_t) instead of 1 / Pi * F_t.
3. Why is the BSSRDF profile multiplied with 2 * Pi * r in ComputeBeamDiffusionBSSRDF? Is this needed for the Catmull-Rom interpolation?
4. Where is the second F_t term in the single scattering integrand? The reference Matlab implementation also seems to lack one of the F_t terms, but it seems to be the other one.
5. Shouldn't cosThetaS (in the single scattering integrand) always be equal to 1? The paper states that theta_i = arccos(-w dot n), so cos(theta_i) = -w dot n. w is the direction of the light beam, which is orthogonal to the surface (and in the opposite direction of the surface normal), so w = -n, so -w dot n = -(-1) = 1.
6. You only seem to consider two intersections when sampling the BSSRDF, but there could be multiple intersections with the surface. It currently is possible that a point sampled using one axis cannot be sampled using a different axis (because it finds a different intersection first), but the PDF seems to assume that the same point is found using all axes. I think you should find all intersections and then randomly pick one, sadly this also requires you to find the intersections for the other axes in order to calculate the PDF.

Cheers, Tom

[tomvbussel]

I took another look at the single scattering, comparing their Matlab code, the paper and the implementation in pbrt-v3. Their Matlab implementation was really confusing as they calculated the same variable in three different ways. In the current implementation you seem to do the same, as cosThetaE and cosThetaS are equal. I simplifided the code a bit:

// Determine angle cosine at exit point
float cosThetaE = t / d;

integral += alpha * PhaseHG(cosThetaE, g) * std::exp(-sig_t * d) *
            cosThetaE / (d * d) *
            (1 - FrDielectric(-cosThetaE, 1, eta));

And with the missing Fresnel term:

// Determine angle cosine at exit point
float cosThetaE = t / d;

integral += alpha * PhaseHG(cosThetaE, g) * std::exp(-sig_t * d) *
            cosThetaE / (d * d) *
            (1 - FrDielectric(-cosThetaE, 1, eta)) *
            (1 - FrDielectric(1, 1, eta));

I rerendered the head scene with the missing Fresnel term included and made a diff of the old and the new image, but the differences are miniscule, but in the ear there seems to be a small difference that cannot be attributed to noise.

I think this solves 4 and 5, but the others still need to be fixed. The authors of "BSSRDF importance sampling" mention 6 in their slides (see slide 49 of these slides).

Cheers, Tom

[wjakob]

Hi Tom,

thanks for your feedback. I am currently doing a pass over the BSSRDF code and will take a look at the remaining points and your suggestion to improve the computation. I'll get back to you in the next few days.

Thanks, Wenzel

[tomvbussel]

Hi Wenzel,

Thanks for the quick reply. I fixed the code in my previous comment, I got confused by the cosThetaE in PhaseHG and the fact that phase function use a different convention (which pbrt-v3 seems to ignore) for the directions than BSDFs.

Cheers, Tom

[tomvbussel]

Hi Wenzel,

I noticed you updated the BSSRDF code (I'm not sure if this is the final version yet), so I decided to have another look at the code. I noticed the new code includes the normalization constant, but sadly it's the wrong constant. The current code evaluates the following BSSRDF:

(1 / Pi^2) * F_t(x_i, w_i) * R(||x_o - x_i||) * F_t(x_o, w_o) / (16 * C_phi(1 / eta)^2)

This is wrong as the normalization constant is squared. To fix this the function S for SeperableBSSRDF should be replaced with:

Spectrum S(const SurfaceInteraction &pi, const Vector3f &wi) {
    return Sw(po.wo) * Sp(pi) * (1 - FrDielectric(CosTheta(wi), 1.f, eta));
}

When sampling the BSSRDF instead of evaluating it, the normalization term and the Fresnel terms seem to be completely ignored. The BSDF (in the case of the KdSubsurface and Subsurface materials) includes one of the Fresnel terms though, so only one Fresnel term is missing.

Personally I would replace the FresnelTransmission in the subsurface materials with a LambertianTransmision (remove the 1 / Pi term from the BSSRDF if you do this) and include both Fresnel terms and the normalization term in Sample_S. Also, when adding the direct subsurface scattering component in the path tracer, the BSSRDF needs to be evaluated again to make sure that the Fresnel term for the right incoming direction is included.

I also figured out why the 2 * Pi * r is (incorrectly) included in the profile in the BSSRDFTable (point 3 in my original post). You need to multiply the profile by this term to compute the diffuse albedo (the constant R_d, see the section "BRDF Approximation" in "A Practical Model for Subsurface Light Transport" by Jensen et al.), but the profile itself should not be multiplied by this term. I think this confusion stems from the fact that R_d is sometimes used to denote the diffuse albedo and sometimes to denote the diffusion profile (in this case it has the scattering point and the light beam direction as parameters).

I still don't understand what all those "sigma_t"s are doing in BSSRDFTable and TabulatedBSSRDF code when sampling the profile (or calculating the pdf of that sample). I assume this has something to do with the fact that you're using the "wrong" sigma_t in ComputeBeamDiffusionBSSRDF, but I'm not sure. I don't fully understand how the table works, as you're using the wrong sigma_tr when tabularizing the profile, the paper seems to suggest a completely different strategy in which they cache tables at shading points.

Finally I'm wondering where you got the sigma_t in the "head" scene from, as it seems to be magnitudes larger than the extinction coefficient for skin reported in "A Practical Model for Subsurface Light Transport" by Jensen et al.

I hope this was helpful.

Cheers, Tom

EDIT: Just ran some tests using the Matlab code that was supplied with the paper. I set g to 0 and eta to 1.33 and sig_a = 0.1, sig_s = 0.9 and sig_a = 0.2, sig_s = 1.8 gave different profiles, despite having the same albedo, so the current tabulation scheme does not work. I think the current tabulation scheme can easily be fixed when using a constant sigma_t, as in this case the sigma_s and sigma_a are uniquely defined for each albedo (given that sigma_t).

[wjakob]

Hi Tom,

thank you again for your comments and suggestions. Here are some responses to your questions:

1. In the initial version, we didn’t account for F_t for the outgoing bounce. Probably not a big deal (this model is extremely approximate in any case). With an architectural change I’m adding right now, it will be easier to incorporate, so the final version will likely have it.
2. Related to (1). Right now, the incident fresnel term is handled by the material’s BSDF, and the outgoing fresnel is neglected (the model just uses a diffuse BRDF lobe). I will see if I can change it, though I don’t think that it will make a significant difference.
3. The 2_pi_r scale factor is needed due to the polar coordinate parameterization of the profile. In this case, one profile entry stores the integral of the profile over the phi parameter. When subsequently integrating over spline segments, one can obtain the integral over the entire BSSRDF profile which provides “effective" material albedo incorporating all orders of scattering. This is an important quantity used in various places (e.g. to solve for the single scattering albedo which causes the material to have a certain effective albedo in the context of texturing) — for this reason, it is is recorded in its own field in the BSSRDFTable.
4. When the BSSRDF sampling routine is invoked, we have already sampled a BRDF that caused the current path to refract into the material. The profile computation assumes that this Fresnel term is already accounted for.
5. I’ve collapsed cosThetaS & cosThetaE into a single constant — good catch!
6. The upcoming version considers all intersections of the probe rays (similar to what is described in the Kulla et al. talk). I don’t believe that that the intersections for the other axes is required for the PDF computation — the technique for each axis (which randomly selects one of the intersections of that axis) is an unbiased estimator on its own. Those techniques can simply be re-weighted using MIS weights based on probabilities for the one-intersection-only case (the resulting estimator will still be unbiased as long as this is consistently done for all techniques).

Your confusion about PhaseHG() is a good indicator that it should probably use the negative sign convention. A corresponding sign flip is now done in HenyeyGreenstein::p() and all other callers of PhaseHG.

About sigma_t: the profile pre-computation is done with sigma_t = 1 (only the albedo varies). When evaluating the profile for a particular sigma_t, this corresponds to a scaling operation that must be applied both to the radius and the actual profile value.

Regarding your last message: I added S() to serve as an explanation in the text (but without being directly called by any part of pbrt). I now plan to return a special instrumented BRDF that calls the BSSRDF profile function. I am still playing with that, so normalizations etc. may be off until this turns into non-dead code.

Cheers, Wenzel

[wjakob]

BTW: can you clarify your experiment? Based on what you said, I would expect that the second profile has the same albedo but is (spatially!) scaled by a factor of 0.5.

[wjakob]

Good catch also about the skin sigma_t coefficients. When texturing was not yet implemented, we rendered the face with Jensen's "Cream" coefficients :) Try "color sigma_t" [772 1050 1490] for "Skin1" instead -- it looks quite a bit nicer too!

The factor of 1000 difference is simply millimeters vs meters.

[tomvbussel]

Hi Wenzel,

Thanks for your answers.

(1 and 2) Okay that makes a lot of sense now. I think that it's trivial though to include this Fresnel term in the single scattering, as it depends on the directions inside the volume, which can be precomputed.

(3) In order to compute the total diffuse reflectance R_d of the BSSRDF you indeed need to integrate the profile over polar coordinates. However when evaluating the profile itself you shouldn't multiply by 2 * Pi * r, which currently is the case.

(6) Currently Pdf_Sp assumes that the same number of intersections is found using all axes, but this doesn't have to be the case. It could very well be the case that when using the surface normal only one intersection is found but when using a tangent multiple intersection are found. In this case Pdf_Sp would overestimate the probability of selecting that point using the tangent. That's the reason why I think that Pdf_Sp should trace a ray for the other axes, to find out the probability of choosing that intersection.

It was clear to me that you're tabulating the profile for sigma_t = 1, but I don't think that's correct. So from what I understand (please correct me if I'm wrong) when looking up a profile, the code first multiplies the distance with sigma_t and then multiplies the result with sigma_t (and clamps that). The problem is that the profile cannot just be spatially scaled, as the profile includes a few terms that non-linearly depend on sigma_tr (which depends on sigma_t).

So regarding my small experiment. I first ran PBD_profile.m using sigma_a = 0.1, sigma_s = 0.9, g = 0, eta = 1.33 and r_array = [1:8]. This gave me the profile for albedo = 0.9 and sigma_t = 1. So let's say I now want to determine the profile for albedo = 0.9 and sigma_t = 2 at r = 1. According to the pbrt code I should first evaluate the tabulated profile at dist = r * sigma_t = 2 and then multiply the resulting value with sigma_t, so the evaluation of the profile with sigma_t = 2 should be 0.0039 * 2 = 0.0078. But when I run PBD_profile.m directly with sigma_a = 0.2 and sigma_s = 1.8 (and the same g, eta and r_array) I get 0.153 at r = 1.

Cheers, Tom

[wjakob]

(3) In order to compute the total diffuse reflectance R_d of the BSSRDF you indeed need to integrate the profile over polar coordinates. However when evaluating the profile itself you shouldn't multiply by 2 * Pi * r, which currently is the case.

Fair enough -- I'll change this in the (unused) definition of S(). However, note that we do want to have this scale factor inside table.profile. This way, we can importance sample an (now 1D-only) distribution over 'r' to get the right spatial distribution. When the scaled profile value is finally divided by the PDF, you can see that the f/p ratio exactly comes out to the effective albedo, which is what we want.

(6) Currently Pdf_Sp assumes that the same number of intersections is found using all axes, but this doesn't have to be the case. It could very well be the case that when using the surface normal only one intersection is found but when using a tangent multiple intersection are found. In this case Pdf_Sp would overestimate the probability of selecting that point using the tangent. That's the reason why I think that Pdf_Sp should trace a ray for the other axes, to find out the probability of choosing that intersection.

Computing the intersections for all axes is really out of the question :). However, keep in mind that this whole projection-based surface parameterization is a heuristic -- there is certainly some leeway in defining a sampling strategy where the current more efficient scheme makes sense, and we should make use of that.

Suppose that there is a strategy (per axis) which always finds a single intersection that will result in a contribution of $f/p * L_i = eff_albedo * L_i $ -- then the current Pdf_Sp() is clearly the right one to use.

This is of course not possible -- in practice, the number $N_k$ of intersections per axis $k$ turns into a random variable. What the current approach is does is integrate this $N_k$ in the final estimator without touching the PDF of the per-axis estimator. This increases the contribution proportional to the number of found intersections.

The problem is that the profile cannot just be spatially scaled, as the profile includes a few terms that non-linearly depend on sigma_tr (which depends on sigma_t).

Respectfully disagreed! It would be very disturbing if the profile depended on sigma_t in a nonlinear way :). Here is more explicit proof: plugging in all definitions and simplifying $\sigma_tr$ yields $\sqrt{3 + \frac{3}{\rho-2}} \sigma_t$ (where $\rho$ is the single scattering albedo)

[tomvbussel]

I'm actually getting:

$\sigma_{tr} = \sqrt{\frac{3\sigma_a(\sigma'_t)^2}{\sigma_a + \sigma'_t}}$

Either way sigma_t does the cancel out in sigma_tr and the profile depends on sigma_tr in a few non-linear ways (a few e^(c * sigma_tr) terms), so the profile depends on sigma_t in a non-linear way. The Matlab code supplied by the authors of the paper seems to agree with me, two profiles with the same albedo, g and eta can't be scaled to fit each other. The paper seems to suggest this as well as they compute a new table at each shading point.

By the way I think the normalization term is still missing from Sample_S.

[wjakob]

You need to set (assuming g=0 to keep things simple) \sigma_s=\rho\sigma_t \sigma_a=(1-\rho)\sigma_t

If you parameterize the profile via the albedo, g and sigma_t, then sigma_t is the scale. If the PBD implementation says differently, then there is a bug in it. There is no nonlinear dependence on scale.

[wjakob]

So I just tried the reference implementation and I see the expected linear behavior (up to minor integration errors)

sigma_a = 0.1; sigma_s = 0.9; g = 0; eta = 1.33; r_array = 0:8;
profile1 = PBD_profile( sigma_a,sigma_s,g,eta,r_array);
spatial_scale = 2;

sigma_a *= spatial_scale;
sigma_s *= spatial_scale;
r_array /= spatial_scale;
profile2 = PBD_profile( sigma_a,sigma_s,g,eta,r_array) / (spatial_scale * spatial_scale);

profile1 ./ profile2

octave-cli-3.8.1:1> test
ans =

   0.99997
   1.00223
   1.02012
   1.03917
   1.00543
   1.00000
   1.00000
   1.00000
   1.00000

This looks correct to me -- maybe you flipped a scale factor in your test?

[wjakob]

By the way I think the normalization term is still missing from Sample_S.

Sample_S() now returns an adaptor BSDF which calls SeparableBSSRDF::Sw(), which has the energy-preserving Fresnel reshaping term (i.e. with the FresnelMoment1() normalization).

[tomvbussel]

Ah, I divided the second profile by the spatial_scale instead of the square of the spatial_scale, as it seems that this is what TabulatedBSSRDF::Sd is doing; first it does a lookup at distance * sigma_t and then it returns fs * sigma_t. Shouldn't the last line of TabulatedBSSRDF::Sd then be:

return (fs * sigma_t * sigma_t).Clamp();

or is something cancelling out during the sampling which makes this unnecessary?

[wjakob]

Hi Tom,

right! While that did not cause issues in practice (the same scale factor occurred in the Pdf_Sd()), it should really have been sigma_t^2 to be consistent with how this function is used.

You were right that including the 2 * Pi * r^2 factor in the profile is problematic. It needs to be in there so that the spline sampling code produces the right distribution in 2D, but Pdf_Sd() and Sd() then need to divide out this term or the weight factors returned by Sample_S() will skew the profile towards higher radii (oh oh..)

I've committed corrections for both issues which will be synced to github shortly. The torso scene visually improved as well -- nice! :)

Thank you for the (very useful!) feedback, Wenzel

[tomvbussel]

Hi Wenzel,

The latest commit looks really good! I think the current code is nearly correct now (I spotted one small bug, see below). The code looks a lot clearer now too. I rendered the head scene again and the visual result looks very pleasing as well. It might be interesting to do another render with a diffuse BSDF (which includes the Fresnel terms), so we can compare the BSSRDF with a BSDF. This should also show if the current implementation is correct.

While that did not cause issues in practice (the same scale factor occurred in the Pdf_Sd()), it should really have been sigma_t^2 to be consistent with how this function is used.

Sd returns a Spectrum, while Pdf_Sd returns a scalar, so these wouldn't exactly cancel out (only when using a single wavelenght sample in the Spectrum). In fact, when rendering with sigma_t instead of sigma_t^2 the skin appears too reddish.

You were right that including the 2 * Pi * r^2 factor in the profile is problematic. It needs to be in there so that the spline sampling code produces the right distribution in 2D, but Pdf_Sd() and Sd() then need to divide out this term or the weight factors returned by Sample_S() will skew the profile towards higher radii (oh oh..)

The current fix is correct mathematically, but sadly it will produce NaNs in Sr and Pdf_Sr at r = 0. These NaNs won't appear in Sample_Sp though, as Sample_Sr won't sample a radius of zero, so you could just keep this bug in there and hope that other readers won't notice. I do agree with you though that this factor needs to be there for the switch to polar coordinates in the sampling (although the 2 * Pi factor actually is produced by the fact that phi is sampled uniformly).

Cheers, Tom

[wjakob]

Hi Tom,

Sd returns a Spectrum, while Pdf_Sd returns a scalar, so these wouldn't exactly cancel out (only when using a single wavelenght sample in the Spectrum). In fact, when rendering with sigma_t instead of sigma_t^2 the skin appears too reddish.

Indeed -- they won't cancel out when sigma_t when varies per channel. I started experimenting and validating the code with just one channel initially, hence this problem was never spotted during development.

The current fix is correct mathematically, but sadly it will produce NaNs in Sr and Pdf_Sr at r = 0. These NaNs won't appear in Sample_Sp though, as Sample_Sr won't sample a radius of zero, so you could just keep this bug in there and hope that other readers won't notice. I do agree with you though that this factor needs to be there for the switch to polar coordinates in the sampling (although the 2 * Pi factor actually is produced by the fact that phi is sampled uniformly).

Hm, that is indeed not so nice (the evaluation routine should clearly never produce NaNs/Infs). I'll make the division conditional on r != 0 and mention in the text that the value is not accurate for r==0.

Assuming you found no other issues, I'll close this ticket. Thank you very much for your feedback -- FWIW you've totally earned a free copy of the book. When the time comes, I'll ping you for your address information.

[tomvbussel]

Hi Wenzel,

Assuming you found no other issues, I'll close this ticket.

Every issue I found was fixed, thanks! Should I find a new bug I will reopen this issue.

Hm, that is indeed not so nice (the evaluation routine should clearly never produce NaNs/Infs). I'll make the division conditional on r != 0 and mention in the text that the value is not accurate for r==0.

Yeah, that should be fine. It's a shame S will return an incorrect zero in that case, but I think that's the best solution, there is no point in completely changing the implementation to fix some dead code that's only there for illustrative purposes.

Thank you very much for your feedback -- FWIW you've totally earned a free copy of the book. When the time comes, I'll ping you for your address information.

Whoa, thanks! I was planning on purchasing the third edition, so that's very nice!

Cheers, Tom

[tomvbussel]

Hi Wenzel,

I noticed another small problem, so I'm reopening this issue. The BSSRDFAdapter currently does not account for the change of medium (from skin to air), while the SpecularTransmission does account for the opposite change of medium (from air to skin), this causes radiance to be incorrectly scaled by a factor of (ext_eta / int_eta)^2. I was also wondering why SpecularTransmission scales radiance with (eta_i / eta_t)^2 instead of the (eta_t / eta_i)^2 which was used in the second edition of the book and was also mentioned by Veach in his thesis.

I also noticed that the second Fresnel transmission term in BeamDiffussionSS is still missing. The paper is a bit vague, but I think there are four Fresnel transmission terms in the case of single scattering: two in S and two in f (which is then multiplied by Q and integrated). Currently the two in S are correctly implemented, but only one of the two Fresnel terms in f (part of the integrand in BeamDiffusionSS) is currently implemented.

Cheers, Tom

[wjakob]

Hi Tom,

the eta^2 factor is actually there: it is SeparableBSSRDF::Sample_Sp(). However, it would be nicer to move it to SeparableBSSRDF::Sw() for symmetry wrt. the dielectric transmission BTDF. (the change will be synced to this repo soon).

Regarding the dielectric: its eta^2 factors look correct to me -- remember that the random walk is done in reverse. When integrating radiance, we're starting at the camera and moving towards the light source -- so you'll need the reciprocal of the scale factor.

I'll look into the BeamDiffussionSS factor later.

Cheers, Wenzel

[tomvbussel]

Oh wow, I didn't notice that the eta^2 factor was already there, my apologies. I would move it to BSSRDFAdapter:f() instead of SeperableBSSRDF::Sw() though, as moving it to Sw() would make SeperableBSSRDF::S() incorrect (unless you fix that too).

I wasn't sure about the factor for transmission, in the second edition of the book (pages 443 and 444) (et_et) / (ei_ei) is used, but I have also seen other sources use (ei_ei)/(et_et). I guess you're right, I just wanted to make sure that it wasn't an accidental typo.

Cheers, Tom

[wjakob]

Good idea, that's even more symmetric. (though I don't think that having the term in S is incorrect, if this is explained in its definition).

[tomvbussel]

It wouldn't be incorrect per se, but it could be confusing as the function S has a very specific meaning in the literature.

[tomvbussel]

Hi,

I was wondering why the current implementation of the BSSRDF ignores bump mapping. The BSSRDF should use the "bumped" positions of the two SurfaceInteractions to preserve the small details that were introduced by the bump map. Currently the Bump function only modifies the surface derivatives, but it should also modify the position of the SurfaceInteraction. I'm not sure if this is intentional, but in my opinion it would have been nicer if the BSSRDF handled it (as the subsurface materials allow it).

Cheers, Tom

[wjakob]

Hi Tom,

of course, it would be fun to add lots of cool extensions :-). However, PBRT is seriously page-bound, meaning: any extra pages we add to describe something must be removed elsewhere. It's a good idea though -- we could make it an exercise at the end of Ch.15.

Cheers, Wenzel

PS: I checked out the Fresnel terms in the SS profile and they look right to me (there is just one transmission term, as the other one is implicitly taken care of by the BSDF).

[tomvbussel]

Hi Wenzel,

That's reasonable. It's just a bit odd that the subsurface materials allow bump mapping while the BSSRDF doesn't support it, but I guess that's completely fine if it's mentioned.

On the Fresnel terms in the SS profile: So the single-scattering portion of the BSSRDF is defined as (you'll have to copy-paste these equations into latex, my apologies):

$S^{(1)}(x_i, \omega_i, x_o, \omega_o) = \frac{1}{\pi} F_t(x_i, \omega_i) R^{(1)}(x_o - x_i) \frac{F_t(x_o, \omegao)}{4C\phi(\frac{1}{\eta})}$

The two Fresnel terms in this equation are currently calculated using the BSDFs. The single-scattering profile is defined as:

$R^{(1)}(x, \omega) = \int_0^\infty r^{(1)}(x, x_r(t)) Q^{(1)}(t) dt$

This is what BeamDiffusionSS estimates. $r^{(1)}$ is defined by equation 13 in the paper. This equation includes another two Fresnel terms. BeamDiffusionSS currently implements one of those two terms, $F_t(\theta_o, \eta)$ if I'm not mistaken, but the other one is missing.

Cheers, Tom

[wjakob]

Hi Tom,

our sampling architecture is slightly different than the PBD paper -- generally the approach is to compute the right physical quantity that makes sense in the context of PBRT rather than meticulously following the definitions in the PBD paper.

We need two fresnel terms for the two refractions at the boundary. I hope you agree.

The first one is taken care of by the FresnelTransmission BxDF.

The second one is more tricky, since the profile tabulates emitted irradiance rather than directionally varying radiance values.

So the approach suggested in PBD is to tabulate the fresnel-scaled irradiance to capture the reduction in exitant radiance due internal reflection at the boundary.

But this only gives a diffuse exitance at surface points, which is a pretty poor approximation -- so PBD adds a third fresnel term to modulate the directional profile. This Fresnel term is normalized and does not reduce the total amount of energy -- it just gets redistributed differently over angles.

Obviously that is a hack, albeit a physically motivated one.

I hope that makes sense :) Wenzel

[tomvbussel]

Hi Wenzel,

My apolgies for my late reply. I was under the assumption that the code was a direct implementation of the paper, so it makes a lot of sense now. Everything looks good now, so I guess I'll close this issue, thanks for your replies.

Cheers, Tom

[Immocat]

Thanks for your @tomvbussel discussion, especially @wjakob's last comment.

Recently(Although in Jan 2018) I've been going into the details about PBD in pbrt v3 book and also the original paper. The third normalized Fresnel term in Sw() has made me confused for two days until I find this page. I really hope that your (@wjakob) last comment was in the book, which would make a lot of sense and save a lot of time :smile_cat:.

Anyway, thanks again for your book and hard work.