path length 1: 3 ways to integrate direct visible light
path length 2: 4 ways to integrate direct lighting
path length 3: 5 ways to integrate 1st bounce lighting
path length 4: 6 ways to integrate 2nd bounce lighting
Simply added all the strategy would result to overbright image (probably 3 times brighter direct visible light, 4 times brighter direct lighting, 5 times brighter 1st bounce lighting....etc) But does simply divide by strategy number work? Something like:
Cs+t=2(ˉx)=1/3(C∗s=0,t=2(ˉx)+C∗s=1,t=1(ˉx)+C∗s=2,t=0(ˉx))
Cs+t=3(ˉx)=1/4(C∗s=0,t=3(ˉx)+C∗s=1,t=2(ˉx)+C∗s=2,t=1(ˉx)+C∗s=3,t=0(ˉx))
Cs+t=4(ˉx)=1/5(C∗s=0,t=4(ˉx)+C∗s=1,t=3(ˉx)+C∗s=2,t=2(ˉx)+C∗s=3,t=1(ˉx)+C∗s=4,t=0(ˉx))
...........
Uhh...... it might work......when you have a camera with intersectable lens + every light in the scene is intersectable + no specular bsdf in the scene, which means: any Dirac distribution light (point/directional/spot lights), bsdf(specular), camera(pinhole camera, orthogonal camera) can break this naive divide method. The reason is, some of the strategy just can't handle certain Dirac distribution: s=0 strategies can't handle point light since ray will never hit the light, t=0 strategies can't handle pinhole camera since ray will never hit the lens, t=1 strategies can't handle direct visible specular surface since fs(ωo,ωi)=0 for every possible ωo,ωi combinations. (We rely on sample bsdf specifically for specular bsdf to get a fs(ωo,ωi) that represent the whole integration result in Dirac distribution case)
Even we assume we won't encounter Dirac distribution, assign equal weight to each strategy doesn't provide the optimal result anyway. We see that certain strategy work particularly good on certain scenario, it doesn't make sense to add in other inferior strategy in this scenario to downgrade the result. The problem is: how we know which strategy work best on which scenario?
Veach95 introduced multiple importance sampling to solve this problem. The sampling combination techniques were applied to several rendering related topics includes: direct lighting sampling, final gathering, and bidirectional path tracing. Lafortune publish bidirectional path tracing around the same time with Veach, but the paper didn't formalize the way to combine different sampling strategies. My personal feeling is that bdpt become a powerful beast after applying multiple importance sampling. It was cool to come up result with multiple strategies, but determine which strategy is better in run time is mind-blowing cool :P
The multiple importance sampling estimator:
F=n∑i=11nini∑j=1ωi(Xi,j)f(Xi,j)pdfi(Xi,j)
where n is the number of sampling techniques, and ni is the number of sample for technique i. With the power heuristic approach Veach propose:
ωi=(nipdfi(Xi,j))2/n∑k=1(nkpdfk(Xk,j))2
For example, we got a length 3 path come up with (s=2,t=2) strategy, and we would like to have its MIS weight. What we need will be figure out the pdf to form this path with other strategies PA(ˉxs=0,t=4), PA(ˉxs=1,t=3), PA(ˉxs=3,t=1), PA(ˉxs=4,t=0) The ni can simply be treated as 1 since we form every strategies for each integration samples. Expand the result:
ωs=2,t=2(ˉx)=PA(ˉxs=2,t=2)2PA(ˉxs=0,t=4)2+PA(ˉxs=1,t=3)2+PA(ˉxs=2,t=2)2+PA(ˉxs=3,t=1)2+PA(ˉxs=4,t=0)2
Veach97 use the notation k=s+t−1 and pi=PA(ˉxs=i,t=s+t−i), rewritten path ˉxs,t as ˉx=x0...xk. This simplify the power heuristic equation to ωs,t=p2s∑ip2s . In above equation:
ωs=2,t=2(ˉx)=p22p20+p21+p22+p23+p24
where:
p0=PA(x3)Pω⊥(x3→x2)G(x2↔x3)Pω⊥(x2→x1)G(x1↔x2)
Pω⊥(x1→x0)G(x0↔x1)
p1=PA(x3)Pω⊥(x3→x2)G(x2↔x3)Pω⊥(x2→x1)G(x1↔x2)PA(x0)
p2=PA(x3)Pω⊥(x3→x2)G(x2↔x3)Pω⊥(x0→x1)G(x0↔x1)PA(x0)
p3=PA(x3)Pω⊥(x1→x2)G(x1↔x2)Pω⊥(x0→x1)G(x0↔x1)PA(x0)
p4=Pω⊥(x2→x3)G(x2↔x3)Pω⊥(x1→x2)G(x1↔x2)
Pω⊥(x0→x1)G(x0↔x1)PA(x0)
All ingredients (Pω⊥, PA, G) can be collected during the path building stage. In fact, Veach even did a clever optimization to avoid duplicate computation between different pi:
ωs,t=p2s∑ip2s=1∑i(pi/ps)2
By treating the current ps=1, we can come up the other pi related to ps using the equation:
p1p0=PA(x0)Pω⊥(x1→x0)G(x1↔x0)
pi+1pi=Pω⊥(xi−1→xi)G(xi−1↔xi)Pω⊥(xi+1→xi)G(xi+1↔xi) for 0<i<k
pk+1pk=Pω⊥(xk−1→xk)G(xk−1↔xk)PA(xk)
and on the reverse side:
p0p1=Pω⊥(x1→x0)G(x1↔x0)PA(x0)
pi−1pi=Pω⊥(xi→xi−1)G(xi−1↔xi)Pω⊥(xi−2→xi−1)G(xi−2↔xi−1) for 1<i≤k
pkpk+1=PA(xk)Pω⊥(xk−1→xk)G(xk−1↔xk)
Pω⊥(x0→x1)G(x0↔x1)PA(x0)
All ingredients (Pω⊥, PA, G) can be collected during the path building stage. In fact, Veach even did a clever optimization to avoid duplicate computation between different pi:
ωs,t=p2s∑ip2s=1∑i(pi/ps)2
By treating the current ps=1, we can come up the other pi related to ps using the equation:
p1p0=PA(x0)Pω⊥(x1→x0)G(x1↔x0)
pi+1pi=Pω⊥(xi−1→xi)G(xi−1↔xi)Pω⊥(xi+1→xi)G(xi+1↔xi) for 0<i<k
pk+1pk=Pω⊥(xk−1→xk)G(xk−1↔xk)PA(xk)
and on the reverse side:
p0p1=Pω⊥(x1→x0)G(x1↔x0)PA(x0)
pi−1pi=Pω⊥(xi→xi−1)G(xi−1↔xi)Pω⊥(xi−2→xi−1)G(xi−2↔xi−1) for 1<i≤k
pkpk+1=PA(xk)Pω⊥(xk−1→xk)G(xk−1↔xk)
There are......some devils hidden in the details that need special treatments.
First, while reusing the path by connecting different path vertex from light path and eye path, the end point pdf need to be reevaluate since the connection changes the outgoing direction:
Second, pk+1pk=0 when pk+1 is a camera with Dirac distribution (pinhole camera) since there are 0 chance you can hit a pinhole camera.
Third, p0p1=0 when p0 is a light with Dirac distribution (point light, spot light....etc) since there are 0 chance you can hit a point/spot light.
Fourth, when xj stands for a surface point with specular bsdf, pj=0 and pj+1=0 Veach97 contains a detail explanation of this exception handling in 10.3.5, the basic idea is, when pj=0:
pjpj−1=Pω⊥(xj−2→xj−1)G(xj−2↔xj−1)Pω⊥(xj→xj−1)G(xj↔xj−1) contains a Pω⊥(xj→xj−1) in denominator which is a number that is not an actual pdf but a hack to get integration working, let's just call it evil alien number for now
pjpj−1pj+1pj=Pω⊥(xj−2→xj−1)G(xj−2↔xj−1)Pω⊥(xj→xj−1)G(xj↔xj−1)Pω⊥(xj−1→xj)G(xj−1↔xj)Pω⊥(xj+1→xj)G(xj↔xj+1) evil alien number still sticks in denominator :|
pjpj−1pj+1pjpj+2pj+1=Pω⊥(xj−2→xj−1)G(xj−2↔xj−1)Pω⊥(xj→xj−1)G(xj↔xj−1)Pω⊥(xj−1→xj)G(xj−1↔xj)Pω⊥(xj+1→xj)G(xj↔xj+1)Pω⊥(xj→xj+1)G(xj↔xj+1)Pω⊥(xj+2→xj+1)G(xj+2↔xj+1)
The interesting thing happen, another evil alien number Pω⊥(xj→xj+1) appears in the nominator. To make thing even more interesting, this two number are actually the same number that can cancel each other. Yes I questioned myself once whether this two number is the same number and the result is true. In pure specular reflection case (mirror) there are both 1. In reflection/transmission combine material like glass, depends on the implementation, it may be 0.5 (giving half half chance to do the transmission or reflection), it may be a F(ωo,ωi) Fresnel number so that you can do a bit importance sampling for this bsdf. 0.5 is self explanatory enough that pdf from ωo to ωi is the same as ωi to ωo. Fresnel number if you examine the equation carefully, pdf from ωo to ωi is still the same as ωi to ωo :)
The conclusion: a specular surface point xj will bring an evil alien number to pdf pj and pj+1 that this two pdf shouldn't be counted in the final weighting calculation, but as long as we move on to pj+2, everything goes back to normal.
With the above gotcha in mind, here is my implementation of multiple importance sampling weight evaluation:
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| float BDPT::evalMIS(const ScenePtr& scene, const CameraPtr& camera, | |
| const std::vector<PathVertex>& lightPath, int s, | |
| const std::vector<PathVertex>& eyePath, int t, | |
| const float GConnect, std::vector<BDPTMISNode>& misNodes) const { | |
| // re-evaluate the pdfForward, pdfBackward for lightPath and eyePath | |
| // end vertex | |
| float pdfSEndForward, pdfSEndBackward; | |
| float pdfTEndForward, pdfTEndBackward; | |
| if (s == 0) { | |
| // eye path end vertex is a light | |
| const Vector3& p = eyePath[t - 1].getPosition(); | |
| const Vector3& n = eyePath[t - 1].getNormal(); | |
| const Light* light = eyePath[t - 1].light; | |
| pdfTEndForward = mPickLightPdf[light->getId()] * | |
| light->pdfPosition(scene, p); | |
| const Vector3& wo = normalize(eyePath[t - 2].getPosition() - p); | |
| pdfTEndBackward = light->pdfDirection(p, n, wo) / dot(n, wo); | |
| // this two are not used in s = 0 case. Set to 0 to be explicit. | |
| pdfSEndForward = 0.0f; | |
| pdfSEndBackward = 0.0f; | |
| } else if (t == 0) { | |
| // light path end vertex is camera lens | |
| const Vector3& p = lightPath[s - 1].getPosition(); | |
| const Vector3& n = lightPath[s - 1].getNormal(); | |
| pdfSEndForward = camera->pdfPosition(p); | |
| const Vector3 wo = normalize(lightPath[s - 2].getPosition() - p); | |
| pdfSEndBackward = camera->pdfDirection(p, wo) / dot(n, wo); | |
| // this two are not used in t == 0 case. Set to 0 to be explicit. | |
| pdfTEndForward = 0.0f; | |
| pdfTEndBackward = 0.0f; | |
| } else { | |
| const PathVertex& sEnd = lightPath[s - 1]; | |
| const PathVertex& tEnd = eyePath[t - 1]; | |
| const Vector3& pSEnd = sEnd.getPosition(); | |
| const Vector3& nSEnd = sEnd.getNormal(); | |
| const Vector3& pTEnd = tEnd.getPosition(); | |
| const Vector3& nTEnd = tEnd.getNormal(); | |
| const Vector3 dSToT = normalize(pTEnd - pSEnd); | |
| if (s == 1) { | |
| float pdfW = sEnd.light->pdfDirection(pSEnd, nSEnd, dSToT); | |
| pdfSEndForward = pdfW / dot(nSEnd, dSToT); | |
| pdfSEndBackward = sEnd.pdfBackward; | |
| } else { | |
| const Vector3 dSEndToSPrev = normalize( | |
| lightPath[s - 2].getPosition() - pSEnd); | |
| pdfSEndForward = sEnd.material->pdf(sEnd.fragment, | |
| dSEndToSPrev, dSToT) / dot(dSToT, nSEnd); | |
| pdfSEndBackward = sEnd.material->pdf(sEnd.fragment, | |
| dSToT, dSEndToSPrev) / dot(dSEndToSPrev, nSEnd); | |
| } | |
| const Vector3 dTToS = -dSToT; | |
| if (t == 1) { | |
| float pdfW = camera->pdfDirection(pTEnd, dTToS); | |
| pdfTEndForward = pdfW / dot(nTEnd, dTToS); | |
| pdfTEndBackward = tEnd.pdfBackward; | |
| } else { | |
| const Vector3 dTEndToTPrev = normalize( | |
| eyePath[t - 2].getPosition() - pTEnd); | |
| pdfTEndForward = tEnd.material->pdf(tEnd.fragment, | |
| dTEndToTPrev, dTToS) / dot(dTToS, nTEnd); | |
| pdfTEndBackward = tEnd.material->pdf(tEnd.fragment, | |
| dTToS, dTEndToTPrev) / dot(dTEndToTPrev, nTEnd); | |
| } | |
| } | |
| // fill the pdf in lightPath and eyePath to misNodes | |
| int k = s + t - 1; | |
| for (int i = 0; i < s - 1; ++i) { | |
| misNodes[i].pTowardLight = i == 0 ? | |
| lightPath[0].pdfBackward : | |
| lightPath[i].pdfBackward * lightPath[i].G; | |
| misNodes[i].pTowardEye = | |
| lightPath[i].pdfForward * lightPath[i + 1].G; | |
| misNodes[i].isSpecular = lightPath[i].isSpecular; | |
| } | |
| if (s > 0) { | |
| misNodes[s - 1].pTowardLight = s == 1 ? | |
| pdfSEndBackward : pdfSEndBackward * lightPath[s - 1].G; | |
| misNodes[s - 1].pTowardEye = ((s - 1) == k) ? | |
| pdfSEndForward : pdfSEndForward * GConnect; | |
| misNodes[s - 1].isSpecular = lightPath[s - 1].isSpecular; | |
| } | |
| for (int i = 0; i < t - 1; ++i) { | |
| misNodes[k - i].pTowardEye = i == 0 ? | |
| eyePath[0].pdfBackward : | |
| eyePath[i].pdfBackward * eyePath[i].G; | |
| misNodes[k - i].pTowardLight = | |
| eyePath[i].pdfForward * eyePath[i + 1].G; | |
| misNodes[k - i].isSpecular = eyePath[i].isSpecular; | |
| } | |
| if (t > 0) { | |
| misNodes[k - (t - 1)].pTowardEye = t == 1 ? | |
| pdfTEndBackward : pdfTEndBackward * eyePath[t - 1].G; | |
| misNodes[k - (t - 1)].pTowardLight = ((t - 1) == k) ? | |
| pdfTEndForward : pdfTEndForward * GConnect; | |
| misNodes[k - (t - 1)].isSpecular = eyePath[t - 1].isSpecular; | |
| } | |
| // iterate from the connection end point to calculate | |
| // relative pdfA and add it to misWeightSum (power heuristic) | |
| float pK = 1.0f; | |
| float misWeightSum = 1.0f; | |
| for (int i = s; i <= k; ++i) { | |
| if (i == 0) { | |
| pK *= misNodes[0].pTowardLight / misNodes[1].pTowardLight; | |
| // exception handling for specular case | |
| if (misNodes[1].isSpecular) { | |
| continue; | |
| } | |
| } else if (i == k) { | |
| // there is no way that light path can hit a delta camera | |
| if (camera->isDelta()) { | |
| break; | |
| } | |
| pK *= misNodes[k - 1].pTowardEye / misNodes[k].pTowardEye; | |
| } else { | |
| pK *= misNodes[i - 1].pTowardEye / misNodes[i + 1].pTowardLight; | |
| // exception handling for specular case | |
| if (misNodes[i].isSpecular || misNodes[i + 1].isSpecular) { | |
| continue; | |
| } | |
| } | |
| misWeightSum += pK * pK; | |
| } | |
| pK = 1.0f; | |
| for (int i = s; i > 0; --i) { | |
| if (i == (k + 1)) { | |
| pK *= misNodes[k].pTowardEye / misNodes[k - 1].pTowardEye; | |
| // exception handling for specular case | |
| if (misNodes[k - 1].isSpecular) { | |
| continue; | |
| } | |
| } else if (i == 1) { | |
| // thre is no way that eye path can hit a delta light | |
| if (lightPath[0].light->isDelta()) { | |
| break; | |
| } | |
| pK *= misNodes[1].pTowardLight / misNodes[0].pTowardLight; | |
| } else { | |
| pK *= misNodes[i].pTowardLight / misNodes[i - 2].pTowardEye; | |
| // exception handling for specular case | |
| if (misNodes[i - 1].isSpecular || misNodes[i - 2].isSpecular) { | |
| continue; | |
| } | |
| } | |
| misWeightSum += pK * pK; | |
| } | |
| return 1.0f / misWeightSum; | |
| } |
And this is how the above first bounce strategies looks like after assign appropriate MIS weight:
You can see the high variance noise (estimator with low pdf) got smooth out and different strategies only preserve the estimator result that with high pdf. The most important thing is: the added on result, still converge correctly compare to the Mitsuba reference.
battle royale round 3: bdpt rendering vs path tracing Mitsuba reference
One thing to notice is that even this scene is not a particularly fun one, bdpt already show that it can use much less samples converging to similar result path tracing did (the above comparison, bdpt use 196 spp, Mitsuba reference uses 400 spp) Though the rendered images are kinda boring (cornell box all the way til the end :P ) The bidirectional path tracing is roughly complete now, and I have to say, OMG!......this journey is definitely not a piece of cake for myself. Time to open up a beer for a small celebration! In the next post, I would try to come up some relatively more interesting images and do a postmortem of this round of implementation.
No comments:
Post a Comment